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CHAPTER  I 


GENERAL  BACKGROUND  AND  SCOPE 


LI  Introductory  overview  of  the  geodetic  boundary  value  problems 

When  studying  the  geodetic  boundary  value  problems  it  is  necessary  to  consider 
theoretical  as  well  as  numerical  aspects  and  limitations.  Historically,  however,  the 
availability  and  nature  of  the  measurements  have  been  the  driving  forces  for  further 
theoretical  considerations  and  development.  It  is  thus  appropriate  not  to  separate  theory 
and  practice,  but  rather,  study  them  as  a  whole  in  a  particular  solution  of  the  geodetic 
boundary  value  problem.  The  evolution  of  geodetic  boundary  value  problems  represents 
the  effort  to  reconcile  theoretical  and  practical  issues,  either  by  making  approximations,  or 
by  varying  the  problem  to  one  with  fewer  approximations,  leading  to  a  progression  of 
problems  of  increasing  complexity. 

The  classical  geodetic  boundary  value  problem  solution  is  the  implementation  of 
fundamental  potential  theory  concepts;  namely  Dirichlet's  principle,  which  states  that 
there  always  exists  a  harmonic  function  (inside  or  outside  a  boundary  surface)  that  takes 
an  arbitrarily  prescribed  set  of  values  on  the  given  boundary.  Furthermore  Stokes' 
theorem  proves  the  uniqueness  of  such  a  harmonic  function  for  a  particular  set  of 
boundary  values.  The  problem  of  determining  a  harmonic  function  outside  a  boundary 
from  its  values  at  the  boundary  is  called  Dirichlet's  problem,  or  first  boundary  value 
problem  of  potential  theory.  Poisson's  integral  is  an  explicit  solution  of  this  problem  for 
a  spherical  boundary  by  means  of  spherical  harmonics  (Heiskanen  and  Moritz,  1967). 

The  third  or  mixed  boundary  value  problem  of  potential  theory  is  the  determination 
of  the  harmonic  function  assuming  boundary  values  which  are  a  linear  combination  of  the 
function  and  its  derivative.  Stokes'  formula  provides  the  solution  to  this  problem  in 
physical  £  desy  where  the  "fundamental  equation  of  physical  geodesy"  is  used  as 
boundar  '  Jition.  It  relates  the  gravity  anomalies  (Ag)  and  the  disturbing  potential  (T) 
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at  the  ellipsoid,  but  its  spherical  approximation  is  used  in  Stokes'  formulation  by 

neglecting  the  flattening  of  the  ellipsoid,  a  sufficiently  accurate  approximation  for  many 
applications.  Still  the  boundary  values  must  be  given  on  the  ellipsoid  which 
approximates  the  geoid.  Since  the  gravity  measurements  are  actually  made  at  the 
topographic  surface  of  the  earth,  they  must  be  appropriately  reduced  to  the  geoid.  This 
becomes  a  serious  shortcoming  of  this  solution  due  to  the  density  assumptions  required  in 
the  various  gravity  reduction  methods,  together  with  the  necessity  of  eliminating  the 
masses  exterior  to  the  ellipsoid,  (Heiskanen  and  Moritz,  1967;  Moritz,  1980). 

In  this  manner  the  anomalous  potential  for  a  "regularized"  earth  is  determined.  The 
determination  of  the  geoid  as  a  level  surface  inside  the  masses,  that  is,  without  a 
regularization  process,  is  made  possible  by  using  two  known  functions  at  the  boundary 
(Molodenskii  et  al.,  1962).  For  example,  Molodenskii's  formula  makes  use  of  the 
gravitational  potential  of  the  external  masses  and  gravity  anomalies  similar  to  free-air 
anomalies.  Still,  the  mass  distribution  assumptions  do  not  allow  for  a  rigorous  solution 
and  the  determination  of  the  geoid  in  this  case  involves  in  addition  an  inverse  problem  of 
potential  theory. 

In  order  to  avoid  the  inverse  problem  and  inherent  mass  distribution  assumptions, 
Molodenskii  reformulated  the  problem  into  a  free  boundary  value  problem;  namely  the 
determination  of  the  earth's  physical  surface  and  the  external  gravitational  field  using  two 
functions  on  the  surface:  the  acceleration  of  gravity  (jf)  and  the  geopotential  (W) 
(Molodenskii  et  al.,  1962).  The  resulting  non-linear  integral  equation  is  linearized  and  is 
expanded  in  powers  of  a  small  parameter  (k),  taking  the  telluroid  as  a  first  approximation 
of  the  unknown  physical  surface  (ibid,  pp.  120-123).  The  solution  of  the  linearized 
integral  equation  is  obtained  through  an  infinite  system  of  integral  equations  which  is 
solved  step  by  step  by  means  of  Stokes'  function.  The  first  equation  gives  a  first 
approximation  which  corresponds  to  Stokes'  formula  for  the  surface  of  the  earth,  while  in 
further  approximations  the  topography  is  considered. 
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Thus,  the  ellipsoidal  heights  of  the  points  on  the  surface  are  computed  from  the  sum 
of  the  height  anomalies  (Q,  determined  from  the  disturbing  potential  on  the  surface,  and 
the  normal  heights  (H),  determined  from  the  values  of  the  geopotential  on  the  surface. 
The  existence  as  well  as  the  uniqueness  of  the  solution  were  examined,  and  although  the 
existence  can  be  guaranteed  for  exact  and  physically  meaningful  data,  the  uniqueness 
cannot  be  guaranteed,  given  gravimetric  data  alone  (ibid.).  Four  additional  parameters 
must  be  determined  by  astronomical  and  geodetic  measurements. 

An  alternative  linearized  solution  to  Molodensldi’s  problem  was  derived  by  Brovar 
(Moritz,  1980)  by  introducing  a  different  harmonic  function  in  the  expression  of  the 
potential  of  a  surface  layer  and  thus  arriving  at  a  different  and  somewhat  simpler  integral 
equation. 


Molodenskii's  problem,  which  is  a  third  type  free  boundary  value  problem  involving 
an  oblique  derivative  has  been  called  the  "modem"  geodetic  boundary  value  problem,  in 
contrast  to  Stokes'  classical  problem.  Certain  simplifications  are  still  present  in  the 
modem  theory.  The  earth  is  considered  rigid  and  uniformly  rotating  around  a  fixed  axis 
passing  through  the  center  of  mass.  Also,  only  the  masses  interior  to  the  earth  are 
considered  and  tidal  effects  are  neglected.  Despite  these  simplifications  the  problem  has 
not  been  solved  in  general  and  recently  there  have  been  important  theoretical  advances  of 
considerable  mathematical  complexity.  In  view  of  the  new  developments,  Molodenskii's 
problem  as  outlined  briefly  above  is  referred  to  as  "the  classical  modem  geodetic 
boundary  value  problem"  (Sansd,  1984). 

A  non-linear  solution  for  Molodenskii's  problem  was  obtained  by  Hormander  using 
Nash's  iteration  method  to  develop  the  inverse  function  theorem  (Moritz,  1980). 

Since  most  of  the  difficulty  is  due  to  the  unknown  boundary,  the  "gravity  space" 
approach  developed  by  Sansd  constitutes  a  significant  contribution,  by  transforming  the 
free  boundary  to  a  fixed  one.  This  is  done  by  describing  the  problem  in  gravity 
coordinates,  that  is  using  the  three  Cartesian  components  of  actual  gravity  =  (gXj,  gX2, 
gx3)  as  the  coordinates  of  the  point  where  they  are  computed.  Since  g  and  W  =  W(]f) 
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are  known  on  the  boundary,  then  the  boundary  is  known  in  the  new  coordinate  system. 
The  serious  disadvantage  of  this  approach  is  the  limitation  to  a  non-rotating  earth,  because 
the  centripetal  acceleration  does  not  allow  for  a  one-to-one  mapping  between  the  ordinary 
Euclidean  space  R3  (x)  and  the  gravity  space  R3  (g).  This  problem  is  handled  by 
considering  only  the  gravitational  potential  and  a  new  function  called  the  adjoint  potential 
is  introduced.  Detailed  reviews  of  the  recent  theoretical  developments  have  been  given  by 
Sansd  (1984)  and  Sacerdote  and  Sansd  (1987). 


A  common  point  between  the  classical  and  the  modem  theory  is  the  requirement  of 
continuous  gravimetric  data  coverage,  although  such  situation  is  far  from  being  realized. 
On  the  other  hand  large  amounts  of  altimetric  data  have  been  available  since  the  last 
decade.  This  prompted  the  formulation  and  study  of  new  geodetic  boundary  value 
problems  of  great  practical  importance,  the  so-called  altimetry-gravimetry  problems. 

A  most  distinguishing  new  element  of  the  altimetry-gravimetry  problems  is  the 
partitioning  of  the  boundary  in  two  parts  where  different  boundary  conditions  hold.  Two 
types  of  altimetry-gravimetry  problems  have  been  proposed.  The  first  one,  presented  by 
Sansd  (1984),  is  a  boundary  value  problem  for  the  Laplace  equation  with  a  Dirichlet 
condition  over  the  sea  surface  (Ss:  fixed)  and  the  geopotential  and  gravity  vector  given 
over  the  land  surface  (Sl:  free).  The  second  one,  presented  by  Holota  (1980),  is  a 
boundary  value  problem  for  the  Laplace  equation  with  a  Neumann  condition  over  the  sea 
surface  (Ss:  fixed)  and  the  geopotential  and  gravity  vector  given  over  the  land  (Sjj  free). 
Only  the  linearized  versions  of  these  problems  have  been  studied  following 
Molodenskii's  classical  treatment,  except  from  adopting  a  different  telluroid  definition  to 
fit  the  known  Ss.  The  relevant  existence  and  uniqueness  theorems  have  been  studied 
under  a  spherical  boundary  approximation  (Sansd,  1983). 

While  the  above  outlined  problems  have  practical  applications,  it  is  currently  the 
case,  and  even  more  so  the  future  situation,  that  different  types  of  data  exist  in 
overlapping  parts  of  the  boundary.  The  corresponding  boundary  value  problems  are 
called  overdetermined.  A  fundamental  distinction  from  the  altimetry-gravimetry 
problems,  where  the  different  observables  cover  disjoint  parts  of  the  boundary  surface,  is 
that  the  redundant  observations  must  be  attached  a  meaning  of  statistical  nature. 
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Obviously,  if  the  observations  were  exact  there  would  be  no  reason  for  retaining 
overlapping  ones.  Hence,  the  solution  of  overdetermined  boundary  value  problems  must 
implement  procedures  capable  of  handling  erroneous  data  in  contrast  to  the  deterministic 
solutions  of  the  classical  and  Molodenskii's  problems.  A  stochastic  approach  to  the 
altimetry-gravimetry  problem  has  been  proposed  by  Bjerhammar  (1983)  utilizing  the 
Gauss-Markov  model  and  the  MINQUE  method  to  estimate  variance  components,  and 
from  them,  the  weight  matrix  to  be  used. 

To  conclude  this  conceptual  overview,  the  principle  of  Integrated  Geodesy  is 
mentioned  as  the  "philosophy"  to  use  all  available  geodetic  (and  other  pertinent) 
observations  in  a  unified  procedure  for  three-dimensional  or  four-dimensional  geodesy 
(Hein,  1986).  In  support  of  the  attainability  of  integrated  data  processing  is  the 
tremendous  advance  in  the  technology  of  the  supercomputers  for  handling  large  amounts 
of  data  and  for  achieving  high  speed  computations. 


1.2  Theory  of  the  overdetermined  geodetic  boundary  value  problems  and 
existing  solutions 


It  is  a  decisive  characteristic  of  the  overdetermined  boundary  value  problem  theory 
that  the  estimation  and  not  the  determination  of  the  gravity  field  is  sought.  A  rigorous 
mathematical  framework  has  been  laid  by  Sacerdote  and  Sansd  (1985)  for  the  solution  of 
overdetermined  boundary  value  problems  in  general,  together  with  the  particular 
formulation  of  geodetically  meaningful  cases,  namely  the  overdetermined  altimetry- 
gravimetry  and  the  gravimetry -gradiometry  problems.  The  essential  ideas  of  this  work 
are  presented  briefly  in  the  following,  starting  from  the  definition  of  the  problem. 


For  the  unknown  function  T,  the  following  conditions  must  be  fulfilled 


V  T=0 

(1)  B1T  =  f,\ 

(2)  B2T=f  2j 


in  Q  :  IR3  =>  Q  < 


on  8Q, 


(1.2-1) 


6 


where  £2  is  the  domain  where  the  Laplacian  (V2)  holds  and  dQ  its  boundary;  that  is,  a 
closed  surface  where  the  boundary  conditions  (1)  and  (2)  are  given.  Bi  and  B2  are 
assumed  to  be  linear  operators,  and  fj  and  f2  the  boundary  data  which  are  functions  > 

defined  on  dQ. 

If  the  problem  defined  by  the  Laplacian  and  boundary  condition  (1)  only  is 
considered  and  assumed  to  have  a  unique  solution  then  the  function  T  is  completely 
determined  and  so  is  f2-  However,  this  is  not  the  case  in  reality  since  fj  and  f2  are 
imperfect  measurements  and  therefore  the  problem  as  defined  has  no  solution.  A  problem 
stated  as  such  is  called  overdetermined  and  a  solution  can  be  obtained  under  a  chosen 
error  minimization  principle.  In  this  sense  a  best  linear  estimate:!  =  (t  1 ,  ^2)  is  found 
from  the  data:  tf°)  =  (fi(0),  f2(0))  by  means  of  the  linear  operator  X,  so  that  the  functions 
fl  and  ^2  are  consistent,  and  a  unique  solution  of  (1.2-1)  is  obtained  from: 

f  =Bt  (1.2-2) 

where  B  a  (Bi,  B2)  and  f  *  X  fl°>. 

Following  these  ideas  one  may  view  the  solution  to  the  overdetermined  boundary 
value  problem  as  the  infinite  dimensional  generalization  of  the  least-squares  estimation  of 
parameters  from  redundant  data. 

Let  yo  be  a  set  of  measurements  of  a  random  variable  with  covariance  matrix  C  and  x 
be  the  parameters  in  the  linear  model  y=Bx.  Then  an  unbiased  estimate  y  is  obtained 
using  the  least-squares  principle: 

(y-yo)TC-1  (9-yo)  =  min  (y  -  yo)T  C1  (y  -  y0),  (1.2-3) 

y=Bx 


where 


9  =  Bx  =  B(BT  C-1B)-1  Bt  C-l  y0. 


(1.2-4) 
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The  equivalent  minimum  principle  for  the  linear  model  (1.2-2)  is  the  infinite  dimensional 
generalization  of  the  principle  (1.2-3),  where  the  vectors  y  and  yo  are  replaced  by  the 
functions  f  and  f(°)  and  the  covariance  matrix  C  becomes  an  operator  C.  Then  the 
variational  principle  may  be  expressed  in  general  as: 

<C't?-f,01,  (f-f10^  =  min  cCff-f10?,  (f-f10^.  (1.2-5) 

f=BT 


The  operator  C  is  symmetric  and  invertible,  defined  from: 
(Cv)i(Xl)  =  JanXCij(Xl>  x2)  vj(X2)dcTx2  •  1*1.2 

j=l 


(1.2-6) 


where  the  kernel  function  Qj  is  the  covariance  function: 

Ci, (X  i,  X 2)  =  E  ([f ,  |X ,)  •  f  j (X  ,|]  [f j (X  j|  -  f j |X 2)]|  (12  7) 

Finally,  the  equation  corresponding  to  equation  (1.2-4)  is  given  as: 

t  =  B(B+C1  B)->  B+C1  fi°>,  (1.2-8) 

where  B+  is  the  adjoint  of  B.  However  this  formula  is  not  accepted  since  the  operator  C*1 
cannot,  in  general,  be  proved  to  be  bounded.  It  is  only  proved  adequate  in  a  simple  case 
where  the  boundary  is  a  sphere  and  the  Dirichlet  and  Neumann  boundary  conditions 
are  given  for  the  data  fi  and  f2  which,  in  this  case,  can  be  expressed  in  spherical 
harmonic  series. 


In  order  to  obtain  an  acceptable  estimator  in  all  cases  Sacerdote  and  Sanso  (ibid.) 
proceed  by  employing  another  minimizing  principle.  The  mean  square  error  principle, 
which  is  proved  equivalent  to  least-squares,  is  written  as: 
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E  (||y  -y||  }  =  min  E^y-yf} 

y=Bx 


(1.2-9) 


where  y  =  E{y).  Imposing  the  condition  of  obtaining  an  unbiased  linear  estimator  the 
above  principle  becomes: 


E  (||y  -  y||2}  =  tr  LTBTQ  1 BLC 


(1.2-10) 


where  L  is  the  unknown  matrix  introduced  in  the  unbiasedness  condition  BLB  =  B,  and 
Q*1  is  any  positive  definite  matrix  used  in  the  norm  definition  in  lRn: 


l|y||2  =  yTQ’1y  (1.2-11) 

The  estimator  finally  obtained  from  the  minimization  of  (1.2-10)  is  independent  of  the 
choice  of  the  norm,  i.e.  the  matrix  Q*1. 

The  above  formulation  is  generalized  to  the  infinite  dimensional  case,  where  the 
result  corresponding  to  (1.2-10)  is: 


E 


£<L+B+Q  ‘BLCe^e^ 

n=l 

trl/B^Q'1  BLC, 


(1.2-12) 


where  (en)  is  an  orthonormal  basis  in  the  Hilbert  space  where  the  norm  is  defined. 


The  convergence  of  the  series  in  equation  (1.2-12)  is  required  and  can  be  secured  by 
choosing  an  appropriate  covariance  operator  C.  In  order  to  assure  a  solution,  the 
minimizing  principle  is  regularized  and  further  modified  to  lead  to  linear  equations. 
However  the  unbiasedness  condition  must  be  released  and  the  new  minimizing  principle 
is: 


,-i 


-i 


t  r  L  B  Q  BLC  +  atr  L  B  Q  BLQ  = 
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min  1 1  r  L+ B  +  Q' 1  BLC  +  atr  L+ B + Q" 1 BLQ 


(1.2-13) 


where  a  is  the  scalar  regularization  parameter.  The  estimator  derived  by  utilizing 
equation  (1.2-13)  is  biased,  but  it  is  proved  that  a  sequence  of  biased  estimates  converges 
to  the  unbiased  one: 


T  =  (b  +  Dq1  b|  B  +  Da1f{0! 


(1.2-14) 


where  Da  =  C  +  aQ  is  proved  in  general  to  be  a  bounded  operator,  except  when  a— >0. 
Comparing  to  equation  (1.2-8)  the  above  result  may  be  viewed  as  a  regularized  least- 
squares  estimator,  since  the  presence  of  the  term  aQ  in  the  definition  of  the  covariance 
operator  is  the  only  difference.  For  the  actual  computation  of  the  estimate  it  is 
suggested  that  a  basis  { nic}  is  introduced  so  that: 


T=XTk"k- 

k=l 


(1.2-15) 


Obviously  the  resulting  system  of  an  infinite  number  of  equations,  involving  the  sum  of 
an  infinite  series,  is  solved  by  truncating  the  series  at  Nmax  and  considering  only  as  many 
equations.  In  conclusion,  the  computational  procedure  is  further  illustrated  in  the 
overdetermined  altimetry-gravimetry  problem: 


BT  = 


®  |  land 
T 

|  sea 


(1.2-16) 


where  the  solid  spherical  harmonics  are  used  as  basis  functions  in  the  T  parameterization. 
Heavy  computation  burden  is  required  for  the  formation  of  the  system  and  also  the 
solution  of  the  system  depending  on  the  choice  of  the  Nmax. 


Because  of  the  extensive  computational  requirements  of  the  global  solution,  Sanso 
(1986)  favors  a  local  solution,  based  on  regionally  defined  functions.  A  basis  of  finite 
element  shape  functions,  for  instance,  would  lead  to  a  patterned  normal  matrix  which 
could  be  exploited  with  sufficient  gain  to  make  the  solution  plausible. 

Other  than  the  computational  difficulties  there  still  remains  a  drawback  in  the  solution 
described  above:  the  dependence  of  the  solution  on  an  arbitrary  regularization  parameter. 
To  avoid  this  Sanso  (1988)  proposed  another  approach,  which  utilizes  a  suitable 
generalization  of  Wiener’s  process  and  Wiener's  integral  over  a  manifold.  On  this  basis, 
the  white  noise  process  is  defined  on  the  boundary,  together  with  its  relation  to  the 
measurement  noise.  Then,  by  defining  the  solution  of  Laplace's  equation  given  white 
noise  as  boundary  values,  a  stochastic  solution  of  a  boundary  value  problem  is  reached 
and  the  minimum  mean  square  error  principle  is  imposed  to  handle  the  overdetermination. 
The  stochastic  model  adopted  in  this  development  allows  for  a  constant  and  independent 
white  noise  in  each  boundary  condition.  As  a  next  step  a  constant  correlation  is  admitted 
between  different  measurements  referring  to  the  same  location.  Further  generalizations 
need  to  be  made  by  considering  variable  variances  and  correlations. 

Despite  the  theoretical  limitations  that  Sans6  has  proved  for  a  least-squares 
adjustment  solution  to  the  overdetermined  boundary  value  problems,  actual  computations 
have  shown  a  rather  successful  practical  result.  In  particular,  high  degree  global 
geopotential  models  (GPM1  and  GPM2)  were  derived  by  Wenzel  (1985),  who  adjusted 
together  mean  gravity  anomalies,  mean  altimetric  geoid  heights  and  satellite  derived 
geopotential  models.  This  procedure  involves  the  solution  of  a  large  system  of  equations 
(of  about  40,000  unknowns)  and  a  large  number  of  observations  (about  100,000).  In 
return,  different  types  of  observations  can  be  included,  without  transforming  from  one 
type  to  another,  while  the  observation  equations  can  be  as  exact  as  possible.  Also,  the 
accuracy  of  the  adjusted  parameters  is  estimated  and  the  residuals  serve  as  good  indicators 
of  the  measurement  errors.  Due  to  the  large  number  of  observations  the  weight  matrix 
(P)  is  assumed  to  be  diagonal  and  the  elements  of  the  normal  matrix  (N)  are  computed  by 
means  of  analytical  expressions,  since  matrix  operations  cannot  be  afforded. 
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Considering  the  estimation  of  the  spherical  harmonic  coefficients  of  a  geopotential 
model  to  degree  and  order  200,  the  normal  matrix  to  be  dealt  with  has  8  x  108  elements  in 
symmetric  storage.  However,  it  is  a  diagonally  dominant  matrix  with  small  off-diagonal 
elements.  In  an  attempt  to  bring  matrix  N  to  a  banded  form,  the  normal  equations  are 
scaled: 

N'  =  SNS  (1.2-17) 

so  that  all  diagonal  elements  of  N'  are  set  to  1.  Then  the  matrix  S  is  a  diagonal  matrix 
such  that: 


[Sii]  = 


1 

yin nr 


(1.2-18) 


All  off-diagonal  elements  [Njj]  <  e,  where  £  is  a  chosen  limit  size,  are  ignored  and  only 
the  ones  [Njj]  £  £  are  taken  into  account.  Since  the  exact  location  of  the  very  small 
[NjjJ's  depends  on  the  data  coverage  and  the  weight  distribution,  it  cannot  be  predicted  a 
priori.  Thus,  to  obtain  a  minimum  bandwidth  of  the  matrix  an  appropriate  ordering  of  the 
unknowns  is  necessary.  As  a  last  step  an  iterative  solution  is  employed  by  using  the 
banded  part  of  the  normal  matrix 


N'  =  B  +  R  (1.2-19) 

where  B  is  the  banded  part  and  R  the  remaining  part  of  N'.  Then  the  iteration  procedure 
is  given  by  the  following  equations: 

li  =  L  -  AXj.i  (1.2-20) 


Xl  =  Xj.i  +  B1  (AT  P  li)  . 


Convergence  is  expected  and  reached  when: 
Xj  -  Xi.i  =  B1  (AT  P  lj)  -»  0  . 


(1.2-21) 
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For  the  computation  of  GPM1  and  GPM2  models,  Wenzel  (ibid)  used  only  the  diagonal 
terms  of  the  normal  matrix  in  B.  Specifically,  for  the  derivation  of  GPM2  complete  to 
degree  200,  55,454  observations  of  mean  gravity  anomalies,  62,003  observations  of 
mean  altimetric  geoid  heights  and  570  potential  coefficients  of  the  GEM  L2  model  were 
used,  a  total  of  1 18,027  observations  and  40,401  unknowns. 

The  problems  that  are  associated  with  this  approach  are  addressed  by  Wenzel  (ibid), 
and  they  are  mainly  the  two  following:  a)  the  spectral  resolution  of  the  data  (especially 
mean  values)  and  b)  the  data  gaps,  which  generate  strong,  short  periodic  variations  in  the 
estimated  model,  unless  a  suitable  stabilization  is  used.  From  the  experience  of  the 
GPM1  and  GPM2  computations,  particular  problem  areas  arc  reported,  which  should  be 
handled  to  avoid  unstable  iterations  in  the  adjustment:  systematic  differences  among 
various  data  sets  should  be  known  k  priori  and  any  existing  data  blunders  should  not  be 
included. 

In  addition  to  these  two  models,  there  exist  a  group  of  high  degree  geopotential 
models  which  are  widely  used  and  can  be  considered  solutions  to  the  overdetermined 
problem,  although  they  do  not  follow  the  typical  definition  presented  here.  On  the  other 
hand,  these  models  are  the  result  of  the  combination  of  all  available  data,  and,  in  that 
sense,  represent  a  non-rigorous  but  practically  accessible  solution.  Among  them  there 
has  been  continuing  refinement  in  the  procedure  as  well  as  the  input  data  sets  used.  For 
example,  the  RAPP81  field,  computed  by  Rapp  (1981)  to  degree  and  order  180  is  based 
on  three  data  sources:  a  1°  x  1°  mean  gravity  anomaly  data  set  derived  from  Seasat 
altimetric  data,  combined  with  a  1°  x  1°  updated  terrestrial  data  set,  yielding  one  set  of 
56,761  values;  also  an  k  priori  satellite  derived  potential  coefficient  set  (GEM  9)  complete 
to  degree  20  with  additional  terms  up  to  (30,  28).  The  combined  anomaly  set  reached 
complete  global  coverage  by  filling  in  the  missing  data  with  values  computed  from  the  a 
priori  coefficient  set.  Overlapping  boundary  values  in  the  two  gravity  anomaly  sets  were 
empirically  merged  to  choose  the  most  accurate  values  and  eliminate  inconsistencies. 
Then,  the  a  priori  coefficients  and  the  global  gravity  anomaly  data  set  were  used  together 
in  a  least-squares  adjustment,  and  the  potential  coefficients  (other  than  the  adjusted  ones) 
were  computed  from  the  adjusted  data  through  an  optimum  quadrature  procedure.  The 
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accuracy  of  the  coefficients  was  estimated  by  error  propagation  considering  measurement 
and  sampling  error. 

As  the  databases  improved  with  time,  more  rigorous  procedures  were  employed  to 
estimate  even  more  detailed  fields.  Rapp  and  Cruz  (1986a)  discuss  the  computation  of 
two  potential  coefficient  sets,  OSU86C  and  OSU86D,  the  only  difference  between  them 
being  that  5,547  geophysically  predicted  gravity  anomalies  were  included  in  the 
derivation  of  the  latter.  The  procedure  followed  was  essentially  similar  to  the  one 
described  above:  a  low  degree  satellite  derived  set  of  coefficients  was  combined  with  a 
global  set  of  1°  x  1°  mean  gravity  anomalies  in  a  least-squares  adjustment.  A  special 
version  of  GEML2  was  used  together  with  six  additional  coefficient  pairs,  after  an 
ellipsoidal  correction  was  applied.  The  anomaly  data  set  was  developed  from  the  merging 
of  the  June  1986  terrestrial  set,  after  applying  downward  continuation,  and  the  Geos 
3/Seasat  altimeter  data  using  empirical  criteria.  The  adjusted  data  were  then  used  in  a 
least-squares  collocation  solution  to  compute  the  remaining  potential  coefficients  complete 
to  degree  250,  together  with  the  associated  accuracy  estimates.  Since  test  runs  of 
replacing  the  collocation  estimation  with  the  integration  formula  procedure  showed  good 
agreement,  the  latter  procedure  was  used  in  deriving  the  next  two  fields,  OSU86E  and 
OSU86F  (Rapp  and  Cruz,  1986  b).  The  adjusted  potential  coefficients  of  the  previous 
fields  were  adopted.  For  the  computation  of  the  remaining  ones  to  degree  and  order  360, 
a  30'  x  30'  gravity  anomaly  data  set  of  149,670  values  was  used  in  an  optimum 
quadrature  procedure.  Accuracy  estimates  could  not  be  easily  obtained  in  this 
methodology  and  it  is  suggested  to  follow  the  values  computed  in  OSU86C/D  models  up 
to  degree  175  and  assign  100%  uncertainty  for  higher  degree  coefficients. 

Comparison  of  these  potential  coefficient  sets  with  GPM1  and  GPM2  and  other 
existing  global  models  with  regard  to  their  theoretical  basis,  problems  and  advantages  is 
given  by  Rapp  (1986).  Further  comparative  analysis  concerning  applications  as  well  as 
the  usefulness  and  future  needs  of  high  degree  models  is  presented  by  Rapp  (1987). 
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1.3  Perspective  and  general  methodology  of  a  collocation  approach 

A  solution  to  the  overdetermined  geodetic  boundary  value  problem,  as  defined  in  the 
previous  section,  is  the  objective  of  this  stuoy.  Considering  the  established  theory,  the 
achieved  practical  solutions  and  the  data  situation,  a  collocation  approach  is  implemented 
in  a  simulation  analysis. 

This  attempt  is  primarily  intended  to  provide  the  methodology  for  the  computation  of 
a  global  gravity  field  model.  The  &  priori  goals  were  the  development  of  a  procedure  as 
general  as  possible  to  handle  all  possible  types  of  existing  information  and  data  types, 
rigorous  with  regard  to  the  solution  of  the  system  and  one-step,  i.e.  avoiding  data  pre¬ 
processing  and  transformation  to  other  types.  Finally  and  most  importantly,  the 
computational  feasibility  of  such  a  procedure  is  sought 

In  general,  an  operational  approach  is  followed  here,  as  opposed  to  a  model 
approach,  eg.  the  classical  boundary  value  problem  solution  (Moritz,  1980,  p.  221).  In 
essence,  operational  is  considered  the  approach  of  solving  a  problem  by  finding  the  best 
possible  solution  in  the  presence  of  all  relevant  information,  instead  of  deriving  a  model 
and  acquiring  the  appropriate  data  to  obtain  its  parameters. 


A  widely  familiar  method  from  "operational  geodesy"  and  "integrated  geodesy",  the 
least-squares  collocation,  is  implemented  here.  Although  collocation  is  mostly  viewed 
with  a  statistical  interpretation,  it  possesses  the  analytical  structure  of  approximation 
methods.  The  particular  properties  of  the  gravity  field  enter  in  a  basic  covariance  function 
from  which  all  other  necessary  covariance  functions  are  derived  via  propagation. 

The  choice  of  the  method  of  collocation  is  advantageous  for  the  following  reasons: 
a)  it  provides  the  best  linear  unbiased  estimate  of  the  predicted  signal  (and  parameters 
when  applicable)  b)  heterogeneous  and  noisy  data  can  be  handled  and  heterogeneous 
signals  can  be  predicted  provided  all  the  necessary  covariance  functions  are  known,  c)  in 
cases  of  data  gaps  the  prediction  reflects  the  corresponding  covariance  functions,  d) 
provides  accuracy  measures  of  the  estimated  quantities. 
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The  most  important  disadvantage  of  the  method  is  the  required  computational  effort, 
which  makes  it  impossible  to  use  in  many  cases.  This  is  because  the  inverse  of  the 
covariance  matrix  of  the  observed  quantities  or  the  solution  of  the  equivalent  system  of 
equations  must  be  computed,  the  former  being  more  time  consuming  than  the  latter. 

Despite  the  large  amounts  of  data  involved  in  the  problem  at  hand,  the 
implementation  of  collocation  is  undertaken  in  this  work  in  light  of  the  following:  a)  the 
block-Toeplitz  structure  of  the  covariance  matrix  of  the  observations  under  certain 
covariance  function  assumptions,  which  is  exploited  in  generating  and  inverting  the 
matrix,  b)  the  sequential  algorithm  of  incorporating  additional  data  groups  based  on  a 
partitioned  inversion,  and  c)  the  recent  advances  in  the  supercomputer  technology  and 
already  achieved  substantial  improvement  over  the  CRAY  X-MP/48  which  was  used  in 
this  work. 


CHAPTER  II 


FORMULATION  OF  A  COLLOCATION  SOLUTION  OF  THE  OVERDETERMINED 

BOUNDARY  VALUE  PROBLEM 


2.1  Fundamental  principles  of  collocation 

The  method  of  collocation  was  introduced  in  geodesy  through  Moritz’  work  of 
gravity  anomaly  interpolation.  Furthermore,  Bjerhammar  presented  the  idea  of 
approximating  the  potential  at  the  points  where  gravity  anomalies  are  measured,  using 
collocation  and  a  set  of  potentials  that  are  regular  down  to  a  sphere  imbedded  within  the 
earth.  It  was  the  valuable  work  of  Krarup  (1969)  which  provided  the  foundation  for  the 
application  of  the  general  collocation  model  in  physical  geodesy.  His  studies  originated 
from  the  instability  suspected  in  Molodenskii’s  boundary  value  problem.  Also,  the  reality 
of  finite  measurements  gave  the  motivation  to  look  at  the  determination  of  the  gravity  field 
as  a  problem  of  interpolation,  or  approximation.  Along  these  lines  it  is  natural  to 
formulate  the  boundary  value  problem  as  an  adjustment  problem,  where  an  improvement 
of  the  boundary  value,  made  minimum  in  some  least-squares  sense,  would  provide  a 
unique  solution.  Krarup  generalized  Moritz’  interpolation  formulation  to  find  directly  the 
potential,  instead  of  using  the  predicted  gravity  anomalies  in  Stokes’  formula.  In 
addition,  his  generalization  included  other  types  of  measurements,  including  satellite- 
related  and  deflections  of  the  vertical,  as  well  as  a  treatment  of  data  error,  in  what  he 
called  a  smoothing  procedure. 


Collocation,  in  applied  mathematics,  is  called  the  determination  of  a  function  by 
fitting  an  analytical  approximation  to  a  set  of  given  linear  functionals  (Moritz,  1980,  p. 
85).  This  definition  is  consistent  with  two  aspects  of  collocation:  the  prediction  aspect 
where  discrete  values  of  the  function  are  predicted,  and  the  one  of  finding  the 
interpolating  continuous  function  as  an  entity.  In  other  words,  there  is  a  finite  and  an 
infinite  dimensional  aspect  of  collocation  (Krarup,  1969,  p.  29).  which  are  both 
applicable  in  physical  geodesy.  The  first  one  when  a  certain  number  of  values  of  linear 
functionals  of  the  gravity  field  are  predicted  at  discrete  points,  and  the  second  one  when 
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the  gravity  potential  is  determined  as  a  function  in  space. 

Consider  a  function  T  in  a  Hilbert  space  H,  and  the  given  linear  functionals  Lk  of  T, 

LkT  =  lk,  k  =  1, 2,  ....  q  (2.1-1) 

where  lk  e  By  defining  the  linear  operator  B  =  [LjJ,  the  above  equation  becomes: 
BT  =  1  (2.1-2) 


with  the  transformation  B:  H  — ►IRQ. 


Assuming  that  a  reproducing  kernel  function  K  (P,  Q)  exists  in  H,  then  it  is  unique 
and  satisfies  two  fundamental  properties: 

K  (P,  Q)  e  H  for  Q  fixed,  (2.1-3) 

f  (Q)  =  <f  (P),  K  (P,  Q)>p  for  all  f  e  H.  (2. 1  -4) 

From  the  above  it  follows  that  the  kernel  function  is  symmetric  and  positive  definite, 
which  is  expressed  by  the  following  relations: 

K(P,Q)  =  K(Q,P)  (2.1-5) 

XXxjXfcKjPi,  PJX), 

i=°k=o  (2.1-6) 

for  all  sets  of  n  points  Pj  and  n  corresponding  numbers  xj  *  0,  where  i  =  1,  2, ... ,  n  for 
all  natural  numbers  n. 

The  kernel  function  may  be  expressed  in  terms  of  the  basis  functions  (pi  (P),  which 
form  a  complete  orthonormal  base  in  H: 


IS 


K(P,Q  =  X  <Pilp)<PilQ)- 

i=l 


(2.1-7) 


It  also  holds  that: 

<pi(P)  =  L?K(P.C8 


(2.1-8) 


and  since  T  e  H,  it  can  be  approximated  by: 


T  (P)  =  X  bi  <Pi  <?)• 

i=i 


(2.1-9) 


where  bj  are  unknown  coefficients,  which  are  determined  by  substituting  (2.1-9)  into 
(2.1-1): 


Lk  T  =  Lk  T  (PJ  =  £  b,  L?K  (pk.  Qi)  =  1k- 

i=l 


(2.1-10) 


Following  the  notation  introduced  in  (2.1-2),  the  base  functions  are  written  as: 


M=bk, 


(2.1-11) 


and  the  equations  (2.1-9)  and  (2. 1- 10)  are  written  as: 


T  =  (BK)  b  and 


(2.1-12) 


B(bk)  b  =  1. 


(2.1-13) 


Solving  (2.1-13)  for  b  and  substituting  into  (2.1-12)  the  function  T  is  determined  from: 
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T  =  (BK)T[B(BK|Y  1.  (2.1-14) 

It  can  be  proved  that  the  solution  given  from  (2.1-14)  minimizes  the  norm  of  T  (i.e. 

1 1  T  1 1 2  ^  II  T  I  1 2)  in  the  Hilbert  space  H.  In  particular  this  norm  has  the  value 
(Moritz,  1980,  p.  214): 

||t||  =<T,T>  =  £X'>ibjL?LjK|P,q. 

1  j  (2.1-15) 


The  above  equation  may  be  expressed  in  the  notation  of  equations  (2.1-1 1)  to  (2.1-14)  as: 


T 


,2 

j  =bTB(BK)Tb, 


(2.1-16) 


and  after  substituting  b  from  (2. 1- 13)  the  norm  is  given  as  function  of  the  data: 


(2.1-17) 


It  is  seen  that  the  determination  of  the  function  T  depends  on  the  choice  of  the  norm  in 
Hilbert  space  H,  and  subsequently  the  kernel  function  K  (P,  Q).  The  particular  choice  of 
the  reproducing  kernel  defines  the  metric  in  the  Hilbert  space,  which  allows  for  a 
geometric  interpretation  of  collocation. 


The  geometry  of  collocation  seen  as  a  least-squares  adjustment  problem  is  shown  by 
(Krarup,  1967,  pp.  34-38),  in  connection  with  the  classical  least-squares  adjustment. 
Two  finite  or  infinite  dimensional  spaces  are  considered  Hi  and  H2.  Given  a  bounded 
operator  A:  Hj  — >  H2  and  an  element  a  e  H2,  find  x  e  Hi:  I  Izl  I2  =  minimum,  where: 

z  =  Ax  -  a  .  (2.1-18) 
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If  the  operator  ATA  is  invertible,  the  unique  solution  of  this  problem  is  the  commonly 
used  least-squares  adjustment  (observation  equations)  formula: 


x  =  (ATA)'1  AT  a . 


(2.1-19) 


Given  a  bounded  operator  B:  H2  — »  Hj  and  an  element  be  Hj,  find  x  e  H2: 


Bx  =  b  (2.1-20) 

and  II  x  1 1 2  =  minimum.  The  solution  in  this  case  is  a  least-squares  collocation  solution: 

x  =  BT  (BBT)-l  b.  (2.1-21) 

The  geometry  of  other  least-squares  adjustment  problems,  a  generalized  adjustment  model 
and  the  collocation  model  is  illustrated  by  Dermanis  (1976).  Most  importantly,  the 
duality  in  interpretation  of  these  models  is  elaborated  and  the  relation  between  the  least- 
squares  adjustment  and  collocation  is  shown. 


In  the  derivations  presented  to  this  point,  the  definition  of  the  metric  in  Hilbert  space 

A 

is  arbitrary,  so  that  there  is  an  infinite  number  of  functions  T  determined  from  equation 
(2.1-14)  which  fit  the  given  data  1,  depending  on  the  choice  of  the  kernel  function  K  (P, 
Q).  In  order  to  obtain  a  unique  and  meaningful  solution,  a  statistical  interpretation  is 
attached  to  the  function  K  (P,  Q):  it  is  identified  with  the  covariance  function  of  T.  This 
is  possible  since  both  functions  share  the  properties  of  symmetry  and  positive 
definiteness.  Also,  all  the  other  relevant  covariance  functions  are  derived  using  the 
propagation  law  of  the  kernel  (or  covariance  propagation  law)  and  the  appropriate  linear 
operators.  Thus,  equations  (2.1-14)  and  (2.1-17)  are  written  as: 


T  —  C’ji  Cjj  1  and 


T— -1 


T  -1  C„  1. 


(2.1-22) 


(2.1-23) 
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where 

C,i  =  cov  (1,1)  and  (2.1-24) 

Cn  =  cov(T,l).  (2.1-25) 


Under  the  statistical  interpretation  this  solution  can  be  derived  by  minimizing  the  variance 
of  the  predicted  quantities  i.e.: 


=  minimum. 


(2.1-26) 


which  corresponds  to  the  minimum  error  norm  approximation: 

|  T  —  T  ||  =  minimum.  (2.1-27) 

An  additional  equation  in  this  case  is  the  covariance  of  the  predicted  quantities: 

-1  T 

Ctt~  CfiCn  C-n  .  (2.1-28) 

The  least-squares  collocation  prediction  method  maintains  the  same  analytical  structure 
and  properties  with  the  analytical  collocation,  while  providing  the  best  linear  estimates,  in 
the  sense  of  minimum  variance  of  the  predicted  signal. 


An  extensive  treatment  of  the  least- squares  collocation  model,  from  the  stochastic 
view  point  is  given  by  Moritz  (1980).  For  the  sake  of  completeness  the  general  least- 
squares  collocation  model  with  parameters  is  briefly  described  here  in  matrix  notation  and 
compared  to  the  observation  equation  model  of  least-squares  adjustment 


Consider  the  linear  model: 

1  =  AX  +  BT  +  n  =  AX  + 1  +  n,  (2.1-29) 
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where  the  data  1  are  represented  by  a  systematic  part  AX,  a  random  signal  part  BT,  and 
random  noise  n.  This  model  reduces  to  least-squares  prediction  with  noise  for  A=0;  for 
B=0,  it  reduces  to  ordinary  least-squares  adjustment.  The  estimation  includes  the  non- 
random  parameters  X  and  the  random  parameters  s  =  [t :  u]T.  The  minimum  principle 
satisfied  is: 

sTC‘^s +  nTCnnn  =  minimum,  (2.1-30) 

and  the  final  estimates  are  given  from: 

^  I  x— i  1  1  x — -t 

x  =  (a1  C  A)  a‘c  1 
s  =  C5tc  '  (l  -  AX) 


(2.1-31) 

(2.1-32) 


where  C  =  Cu  +  C„n-  1°  contrast,  the  minimum  principle  involved  in  least- squares 
adjustment  is: 


T  -1 

n  Cnnn  =  minimum, 


(2.1-33) 


where  n  is  the  only  random  quantity.  In  the  case  of  prediction  and  filtering  (or 
smoothing),  equation  (2.1-32)  becomes: 

-1 

s  =  Cst(C,t+Cnn)  1.  (2.1-34) 


The  above  estimates  are  optimal  in  the  minimum  variance  sense,  provided  the  covariances 
involved  are  known.  When  0^=0,  equation  (2.1-34)  reduces  to  equation  §  =  Cst  C"u  1, 
known  as  Wiener-Kolmogorov  formula  (ibid.,  p.  211). 
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Other  than  the  derivations  outlined  above,  these  formulae  can  be  obtained  in  a 
deterministic  way,  by  adopting  quadratic  norms  and  solving  the  corresponding  variational 
problems. 

2.2  Estimation  of  spherical  harmonic  coefficients  for  the  disturbing 
potential. 

The  aim  of  this  work  being  to  find  the  function  of  the  anomalous  gravity  potential  T 
of  the  earth,  a  set  of  basis  functions  is  needed  so  T  can  be  expressed  as  their  linear 
combination.  Then  the  unknown  coefficients  are  to  be  calculated  as  shown  in  equations 
(2.1-9)  and  (2.1-10).  The  foundation  study  of  Krarup  (1969)  has  discussed  the 
usefulness  of  spherical  harmonics  as  such  a  base  from  the  practical  point  of  view,  beyond 
the  fact  that  they  represent  a  solution  to  Laplace’s  equation.  These  concepts  will  be 
briefly  exposed  in  the  following. 

Let  £  be  the  space  outside  a  sphere  with  radius  R  and  surface  o  (a  e  £)  and  sets  of 
potentials  <p  which  are  regular  in  £  and  at  infinity  so  that: 

lim  (p  (P)  =  0.  (2.2-1) 

P-4  OO 


A  set  S  of  these  potentials  which  are  continuous  in  £  and  o  have  continuous 
boundary  values  on  the  surface  a  of  the  sphere.  It  is  proved  that  if  the  boundary  values 
are  known,  then  <p  can  be  found  by  means  of  Poisson's  integral  and  that  Poisson's  kernel 
can  be  expressed  in  spherical  harmonics.  From  that,  a  symmetric  kernel  function  K  (P, 
Q)  can  be  defined  as: 


K(P,Q)=£(2n+l| 

n=0 


.2  \ 


n+1 


rPrQ  / 


Pn(cos  y) 


(2.2-2) 
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It  is  shown  (ibid.  p.  44)  that  K  (P,  Q)  is  the  reproducing  kernel  for  a  Hilbert  space  H, 
which  can  be  proved  to  consist  of  potentials  regular  in  £,  with  square  integrable 
boundary  values  on  the  sphere.  After  defining  the  functions  {<pn  (P)}  by: 


<P„(P)  = 


1,0 

n+1  \ 

M 

|  ^nm  (Qp>  ^>p)>  y 

\  /  V 

n+1  / 

Snm(Qp>  ^-p)>  m<0  1 

(2.2-3) 


where  R.nm  and  Snm  are  fully  normalized  surface  spherical  harmonics,  then  (2.2-2)  is 
written: 


i*"t  J*,  m  m 

k(p.q=X  X  <p„  <p„  tot, 

n=0  m=-n 


(2.2-4) 


which  means  that  the  spherical  harmonics  {tp™  )  form  a  complete  orthonormal  system  for 
the  Hilbert  space  H.  Then  for  every  potential  <p  e  H  a  series  expansion  may  be  written: 


”  "  mm 

<P  IP)  =  X  2,  an  <P„  (P).  for  Pe  I. 

n=0  m=-n 


(2.2-5) 


There  exists  a  sphere  of  radius  Ri,  for  which  the  above  series  converges  uniformly  on  the 
surface  and  outside  of  any  concentric  sphere  R,  such  that  R>Ri.  Additionally,  there  are 
regions  of  convergence  within  the  limit  sphere  Rj.  In  this  sense  the  series  provides  an 
analytical  continuation  of  the  exterior  potential  within  the  body.  The  potential  of  a 
homogeneous  sphere  for  example,  can  be  analytically  continued  into  the  whole  space 
except  the  center.  In  the  current  application  it  is  of  interest  to  know  whether  the  limit 
sphere  for  the  disturbing  potential  would  be  located  below  the  earth’s  surface.  It  can  be 
proved  that  even  if  that  were  true,  the  addition  of  some  mass  at  a  location  above  the 
surface  will  result  in  a  potential  with  singularity  at  this  particular  location.  Thus,  it  is  not 
practically  meaningful  to  speak  of  convergence  of  the  series  near  the  surface,  since  it  is  a 
very  unstable  condition.  As  it  is  stated  in  the  literature  (ibid.  p.  51,  Moritz,  1980,  p.  64): 
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"if  the  series  were  convergent  at  the  surface  of  the  earth,  a  displacement  of  a  single  grain 
of  sand  would  spoil  the  convergence".  The  remedy  to  this  problem  is  supplied  by 
Runge's  theorem,  which  states  (Krarup,  1969,  p.  54):  "Given  any  potential  regular 
outside  the  surface  of  the  earth  and  any  sphere  in  the  interior  of  the  earth.  For  every 
closed  surface  surrounding  the  earth,  there  exists  a  sequence  of  potentials  regular  in  the 
whole  space  outside  the  given  sphere  and  uniformly  converging  to  the  given  potential  on 
and  outside  the  given  surface".  Hence,  the  potential  can  be  approximated  arbitrarily  well 
by  polynomials  of  spherical  harmonics  down  to  the  surface  of  a  sphere  interior  to  the 
earth. 

In  order  to  proceed  with  finding  the  spherical  harmonic  coefficients  of  the 
approximation  of  the  disturbing  potential  using  collocation,  a  choice  of  the  reproducing 
kernel  function  K  (P,  Q)  is  necessary  to  obtain  a  unique  solution.  This  is  accomplished 
by  introducing  a  stochastic  interpretation  of  the  kernel  as  the  covariance  function  of  T. 
This  is,  however,  faced  with  both  conceptual  and  practical  problems. 

First,  the  covariance  function  cannot  be  found  from  formal  stochastic  theory  directly, 
since  this  would  require  a  population  of  earths  and  the  corresponding  number  of 
realizations  of  T  measured  at  all  points.  Considering  T  as  a  stochastic  process  on  the 
sphere,  an  empirical  covariance  function  may  be  computed  from  only  one  realization  of 
the  process  by  replacing  the  expectation  operator  E  (i.e.  phase  average)  by  a  suitable 
space  average  M,  and  assuming  that  the  process  is  an  ergodic  process  on  the  sphere 
(Moritz,  1980,  p.  285). 

A  second  problem  is  that  T  is  not  directly  measured.  This  is  handled  with  the 
propagation  law  of  covariances  which  allows  for  the  covariance  K  (P,  Q)  to  be  computed 
from  the  covariances  of  its  functionals  (e.g.  C  (P,  Q)  of  gravity  anomalies). 


Despite  these  problems,  a  positive  consideration  is  that  certain  properties  of  K  (P,  Q) 
are  known,  namely  its  symmetry,  positive  definiteness  and  harmonicity,  as  discussed  in 
Section  2.1. 
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Finally,  another  assumption  is  made  considering  T  to  be  a  rotation  invariant 
stochastic  process.  This  results  in  a  covariance  function  which  depends  on  the  spherical 
distance  between  the  points  P  and  Q.  Such  a  function  is  called  homogeneous  (origin 
independent)  and  isotropic  (azimuth  independent)  (ibid.  p.  283).  Since  the  earth  is  not  a 
sphere,  such  assumption  may  be  questionable.  It  is  nevertheless  considered  a  feasible 
approximation  of  the  statistical  model  (Krarup,  1969,  p.  22). 


In  the  following  the  formulae  for  the  least-squares  collocation  prediction  of  the 
spherical  harmonic  coefficients  are  presented.  The  procedure  is  essentially  the  one 
outlined  by  Moritz  (1980,  sec.  21).  Certain  other  equations  are  derived  by  Sjoberg 
(1978). 


By  eliminating  the  systematic  part  from  equation  (2. 1  -29)  the  measurements  are 
modelled  as: 


1  =  t  +  n,  (2.2-6) 

where  1  are  the  measurements,  consisting,  for  this  study,  of  gravity  anomaly  (Ag)  and 
undulation  (N)  data.  The  signal  part  of  1  is  t  and  n  is  the  measurement  noise.  Since  the 
disturbing  potential  function  T  e  H  is  expanded  into  an  infinite  series,  the  spherical 
harmonic  coefficients  comprise  an  infinite  dimensional  vector, 

s  =  [sj  S2  S3  ....  ]T  ,  (2.2-7) 

of  the  Hilbert  space  I2.  The  equivalence  between  T  and  s  defines  an  isomorphism 
(Moritz,  1980,  p.  213)  between  the  corresponding  spaces: 

T  -*  s  :  H-M2,  (2.2-8) 

which  is  also  shown  to  be  isometric: 

Ihf-w2- 


(2.2-9) 
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Thus,  equation  (2.2-6)  is  written: 


1  =  Bs  +  n. 


(2.2-10) 


where  the  operator  B  is  expressed  as  a  q  x  «  matrix  (where  q  is  the  dimension  of  the 
vector  1).  The  kernel  K  is  a  symmetric  infinite  matrix,  and  the  formulae  of  the  finite 
dimensional  case  may  be  applied: 


s  -  Cst(Ctt  +  Cnn)  1, 


(2.2-11) 


under  the  condition: 


sT  KSs  +  nT  D1  n  =  minimum 


(2.2-12) 


where: 


Cst=  cov  (s,  t)  =  KB  * 

Ctt=  cov(t,  t)  =  BKB1 
Cnn  =  cov  (n,  n)  =  D 
CSJ  =  cov  (s,  s)  =  K 


(2.2-13) 


The  accuracy  estimate  of  the  predicted  signal  is  given  by: 

Ess  =  Qs  "  Cst  (C„  +  Cnn)‘^  Cst^  •  (2.2-14) 

Next,  the  formulae  for  the  computation  of  the  covariances  in  (2.2-13)  are  given,  starting 
from  the  basic  covariance  function  K  (P,  Q).  A  detailed  derivation  is  given  by  Rapp 
(1988). 


K  (P,  Q)  =  M  {T  (P)  T(Q)}  =  K  (Vpq) 
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—  flK  r  f  2K  t(0,x)  T  (0',X')  sin  0d0dXdcc 
gK2  •'juo  *'e=o*'a=o 


(2.2-15) 


where 


cos  y  =  cos  0  cos  0'  +  sin  0  sin  0'  cos 


(v-x) 


(2.2-16) 


(2.2-17) 


In  the  above  expression  Cnm  and  Snm  are  the  conventional  fully  normalized  spherical 
harmonic  coefficients  of  the  disturbing  potential,  where  R  is  the  associated  scaling 
parameter.  Also  Pnm  are  the  fully  normalized  associated  Legendre  functions. 


After  performing  the  integration  in  equation  (2.2-15)  (Heiskanen  and  Moritz,  1967,  pp. 
257-258)  and  considering  no  zero  and  first  degree  harmonics: 


K|P,Q)=Xk„P„(cos¥) 

n=2 


(2.2-18) 


where  kn  are  the  potential  degree  variances  given  by 


(2.2-19) 


Cnm  and  Snm  are  fully  normalized  spherical  harmonic  coefficients.  The  extension  of 
K(P,Q)  in  three  dimensions  is  done  through: 


k(p.qi=IM 

n=2 


.2  \ 


rplQ/ 


n+1 


Pn(cos  y). 


(2.2-20) 
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Now,  the  potential  is  expressed  in  terms  of  Laplace  surface  harmonics  Tn  (0,  X): 

~  ,_,n+l 

T(P|=Z|rj  T„(e,x). 

n=2\  p  (2.2-21) 

and  the  radial  component  of  the  gravity  anomaly  at  a  point  is  also  expressed  in  terms  of 
Tn  (6,  X)  as: 


Ag(P)  =  l 

p 


°°  /p,n+l 

5>-l|£  T„(8,X). 

n=2  V rP/ 


(2.2-22) 


Proceeding  with  the  same  argument  as  the  one  used  for  the  derivation  of  K  (P,  Q),  the 
covariance  function  for  the  gravity  anomalies  is  obtained: 


/  r2  \ 


cov(Agp,  AgJ=C(P,Q)=  £cn  — —  Pn(cosy), 


n=2 


n+2 


\rPrQ  I 


(2.2-23) 


where  cn  are  the  gravity  anomaly  degree  variances  given  by  (Rapp,  1988,  p.18) 


c 


n 


(2.2-24) 


By  substituting  equation  (2.2-19)  in  the  above  and  making  a  spherical  approximation  for 
y  evaluated  at  a  sphere  of  radius  R,  (i.e.  y=  kM/R2)  the  above  equation  becomes: 


c„  =  yVi)2£(cL*sL). 


m=0 


(2.2-25) 


Note  that  the  anomaly  degree  variances  evaluated  from  the  above  equation  refer  to  a 
sphere  of  radius  R  (Rapp,  1988,  eq.  (101)).  If  used  for  the  computation  of  covariance 
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between  data  points  on  the  surface  of  a  sphere  with  radius  Re,  the  equation  (2.2-23)  is 
usually  written  as: 


C(y)=  £  cnsn+2pn(cos  y), 

n=2  (2.2-26) 


where 


s  = 


The  same  result  for  C  (\j /)  is  reached  by  using  the  covariance  propagation  formula,  as  is 
illustrated  in  the  derivation  of  the  covariance  function  between  the  anomaly  and 
undulation  signals.  In  general  it  holds: 


covj/lgp,  N(J=L,PLjQK|P,Q). 


(2.2-27) 


Since,  without  loss  of  generality  the  radial  component  of  the  gravity  anomaly  can  be 
expressed  in  a  spherical  approximation  as: 


a  dT 


-  — T 
r„  p 


the  operator  L*  is  defined  as 


Lf  = 


a__2 
ar  r 


(2.2-28) 


and  from  Nq  =  To/y,  the  other  operator  involved  is 


lq-I 

1  Y 


(2.2-29) 
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By  applying  the  above  operators  to  K(P,  Q),  using  a  constant  y.  and  further  manipulation: 


cov 


(AgpN^i 

y 


SKjP.Q) 

-1k(p,qi 

dr 

p  rp  J 

n+1 


Pn(cos  v). 


(2.2-30) 


For  computations  of  covariance  between  points  on  the  sphere  of  radius  Re,  the  above 
formula  yields: 


covjAgp,  nJ  =  —  £-ysn+1  pn(cos  v)  • 
yREn=2n 


(2.2-31) 


Similarly,  the  autocovariance  function  for  the  undulation  measurements  is  obtained  as: 


“  /  1  \"+1 

“vINp.N^-irlkJ-S-  P„(cosV), 
„n-2  \TPTQ> 


(2.2-32) 


and  for  computations  on  the  sphere  Re: 


COvfNp,  ^TS^PnfcOS  v) 


n=2 (n-1) 


(2.2-33) 


When  including  data  which  are  area  means  of  point  data  values  the  need  for  mean 
covariance  functions  arises.  For  a  less  time-consuming  computation  of  such  covariances 
for  signals  on  the  sphere,  the  equations  (2.2-26),  (2.2-31)  and  (2.2-33)  are  modified  by 
introducing  a  smoothing  operator,  called  the  Pellinen  operator  (Sjoberg,  1978)  (in-  The 
corresponding  equations  are: 


(2.2-34) 


C(P.Q|=£m?c11s"2P11(cosv), 

n=2 


COvjAgp,  N(J=-—  X  PnPn~Tsn+lpn(cosvl 

yR  £  n=2  n  1 


COv(Np,  Nq)  =  — XPnPn— -^sn+lpn(c°S  v). 

y 


^  •  ii  •  it  2 

-2  M 


where. 


B  =_1 _ i_ 

"  1-COS  Vo  2n+1 


Pn-1  (cOS  Vo)  *  Pn+ 1  (cOS  V0) 


(2.2-35) 


(2.2-36) 


and  VO  is  the  equivalent  spherical  cap  radius  for  the  block.  When  equal  area  means  are 
used,  as  given  by  Sjoberg  (1978),  the  value  of  VO  is  constant,  and  for  blocks  at  the 
equator  is  given  from  the  equation: 


sin 


l*i  ,8*AX 
471  / 


(2.2-37) 


and  in  equations  (2.2-34  to  36)  it  holds  pn  =  P  n  •  For  area  means  referring  to  an 
equiangular  grid  the  value  of  VO  is  computed  as  shown  by  Katsambalos  (1979), 


cos  Vq=  1  “ 


8  sin 


<pN-sin  9S) 


2ti 


(2.2-38) 


For  efficient  computation  of  the  Pn  factors  the  following  recursive  formula  is  used: 
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2n-l  n-2 

P»=-s?rcosvo|3»-'-iiTrp"-2’ 


where  the  starter  values  are  Po  =  1  and  Pi  =  (1  +  cos  yo)/2- 


(2.2-39) 


A  second  group  of  required  covariances  pertains  to  the  cross-covariances  (Cst) 
between  the  observed  signal  and  the  predicted  spherical  harmonic  coefficients,  defined  by 
the  relations  (2.2-13).  The  matrix  K,  which  defines  the  metric  in  the  space  I2,  may  be 
derived  from  the  function  K  (P,  Q)  by  propagation.  Starting  from  the  representation  of 
the  anomalous  potential  on  the  sphere  Re: 


°°  n  /  D  1  1 

r(R&e.x)=  X  £  hr  s™Rnm(8A)+b„msnn,(e,x)L 

n=2m=0'  E' 


(2.2-40) 


the  coefficients  are  found  from. 


anm=—  f  t(8,A.)r nm(0,x)da,  and 
4  k  J, 

=  —  f  T(e,x)s„m(e,x)da 

4tt  J  . 


(2.2-41) 


Applying  covariance  propagation  to  K  (P,  Q)  and  using  the  operators  defined  in  (2.2-41) 


- l—j  j  K|P.Q|Rm,(e,X)Rpq(0,X)dad<5' 


o  "  O' 


(2.2-42) 


where  do  and  do7  are  differential  surface  elements  on  the  sphere.  Next,  the  function  K 
(P,  Q)  is  substituted  from  equation  (2.2-17)  where  Pn  (cos  \jr)  is  expressed  by  the 
decomposition  formula.  After  the  interchange  of  summation  and  integration  and  applying 
the  orthogonality  properties,  the  result  is  (Moritz,  1980,  p.  160) 
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cov  (ann,  anmj  =  — — -  and 


2n+l 


cov  (anm,  apj  =  0,  for  n*p  or  m*q  or  both. 


Similarly: 


covjbrmr  bnm)  = 

cov  (bnm>  bpq)  =  0,  for  n*p  or  m*q 
cov  (innv  bpj  =  0,  always. 


or  both. 


(2.2-43) 


(2.2-44) 

(2.2-45) 


For  the  usual  representation  of  T, 


T(V,e,X)  =  ^£(y  X  [cnmc°smX  + Snmsin  mX  Pnm(cos6), 

n=2 '  *  '  m=0 


(2.2-46) 


the  corresponding  covariance  among  the  coefficients  Cnm  and  Snm  is  easily  found  as: 


2 

cov  (Cnm.  C^)  =  COV  fSnm>  S^)  =  2n+\  ’ 


(2.2-47) 


and  zero  in  any  other  coefficient  combination.  Then  the  matrix  K  is  a  diagonal  matrix  of 
infinite  dimension: 


C„=K  = 


:n  0 


k22 


0 


.33 


(2.2-48) 


where  the  value  of  kji  is  the  same  for  all  coefficients  of  a  particular  degree  n  and 
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k„„  =  ■ 


Y  (n-l)  (2n+l) 


=w- 


(2.2-49) 


In  order  to  derive  the  matrix  B  defined  in  (2.2-10),  the  gravity  anomaly  and  the 
undulation  are  written  in  terms  of  spherical  harmonics: 

~  ,„,n+2  n  r 

Ag(r,0,X)  =  y£(n-l)  -  £  Cnmcos  mX  +  S^sin  mX  P^fcosO), 

n=2  '  ’  m=0 


(2.2-50) 


and 


N(r,8A)  =  Rlif.  2 


n=2  '  '  m=0 


Cnmcos  mX  +  S^sin  mX 


Pnm(cose). 


(2.2-51) 


Hence  the  element  of  the  B  matrix  relating  the  spherical  harmonic  coefficients  of  degree 
and  order  (n,  m)  to  the  anomaly  measured  at  point  P  is: 


B. 


,n+2 

=  Y(n-l)j-J  cosmXpP^fcoseJ 


(2.2-52) 


B. 


(P  ,n+2 

-j  sinmXpPnm(cos0p). 


(2.2-53) 


Equations  corresponding  to  (2.2-50)  and  (2.2-51)  for  mean  data  are: 


_  v  **  /t?  \n+^  n  "  r 

Ag  (r,0,X)  =  —  £  (n-l)l—  £  C nmJ  P^jcos  0)  cos  mX  da 

Cfn=2  ' 1  m=oL  ° 
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+Snm  Pnm(cos  e)  sin  mX  daj 

cnm  /  pnm(cos  0)cos  mX  da  + 
•'o 

+  S  nm  J  Pnm  (cos  e )  sin  mX  da  , 

•'a 


(2.2-54) 


(2.2-55) 


where  a  is  the  area  of  the  block  on  the  unit  sphere,  with  block  boundaries  the  parallels 
<Pn,  <Ps  and  the  meridians  Xe,  X\v,  given  by: 


a  =  (XE-  Xw)  (sin  cpN  -  sin  <ps) . 


(2.2-56) 


Note  that  R  in  the  above  equations  is  the  scaling  parameter  associated  with  the  coefficients 
Cnm  and  Snm  .  The  integrals  are  computed  analytically  on  the  unit  sphere,  where  da  = 
sin9  d8  dX,  and  using  the  notation  Plnm  for  the  integrated  Legendre  functions: 


PInm=  f  pnm(cos8)sin0de, 
*/0=0, 


(2.2-57) 


the  integrals  are  as  follows: 


—  PInmr  I 

integral  term  of  Cnm  =  ~~~  sin  mXE  -  sin  mX^J,  m?*0 


integral  term  of  Cn0  =  P^nm[^E-  ^w) 


integral  term  of  S  nm  = 


PI, 


m 


cos  mXw-  cos  mXE,  m *0 


integral  term  of  Sno  =  0 . 


(2.2-58) 

(2.2-59) 


(2.2-60) 

(2.2-61) 
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After  further  manipulation  the  above  equations  give  the  elements  of  the  Cst  matrix,  which 
are  also  derived  by  Sjoberg  (1978,  p.  8),  although  following  a  different  argument. 


A  comparison  between  the  practical  application  of  the  collocation  formulation 
presented  here  and  the  least-squares  adjustment  is  made  by  Moritz  (1980,  pp.  166-167). 
By  estimating  a  finite  number  of  coefficients  s  =  [si  S2  ...  sn]t,  the  infinite  matrix  B 
reduces  to  dimension  q  x  N;  for  N<q  an  overdetermined  system  of  equations  is  to  be 
solved.  The  solution  equation  (2.2-1 1)  is  written: 


bkbt+d 


(2.2-62) 


and,  by  applying  a  matrix  identity: 


T  -1 

B  D  1. 


(2.2-63) 


Since  the  matrices  K  and  D  are  diagonal,  the  difference  between  the  implementation  of 
equations  (2.2-62)  and  (2.2-63)  is  that  the  first  requires  the  inversion  of  a  q  x  q  matrix, 
whereas  the  second  of  an  N  x  N  matrix.  If  the  matrix  K  has  the  form: 


K  =  X\ 


(2.2-64) 


where  I  is  the  unit  matrix  and  X.  a  scalar,  the  equation  (2.2-63)  becomes: 


s  =/bTD  1  B  +  —  i)  BTD'll. 


(2.2-65) 


and  for  X,  -» 


s  =  (btd  1  b)  BTD'1  I. 


(2.2-66) 
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which  is  the  least-squares  adjustment  solution  of  the  observation  equation  model  (2.2-10) 
under  the  condition  nT  D_1  n  =  minimum  and  treating  the  coefficients  s  as  non-random 
parameters. 


2.3  Toeplitz  pattern  and  inversion  algorithms 

Practical  applications  of  collocation  are  often  limited  due  to  the  tedious  computations 
required.  Most  of  the  effort  is  spent  in  formulating  and  inverting  the  covariance  matrix  of 
the  observations,  which  makes  the  method  prohibitive  to  use  in  applications  where  a  very 
large  amount  of  measurements  is  involved.  In  order  to  work  around  this  problem, 
studies  have  been  made  in  various  cases  to  identify  special  properties  of  the  covariance 
functions  and  the  underlying  process,  that  result  in  patterned  matrices  consuming  less 
computational  time  for  their  formation  or  inversion  or  both. 

In  geodetic  applications  a  certain  pattern  develops  from  data  regularly  sampled  on  a 
sphere  (Colombo,  1979).  For  data  given  on  an  equiangular  grid  under  the  assumptions 
of  complete  coverage  of  the  sphere,  excluding  the  poles,  and  the  covariance  func'ion 
dependent  only  on  the  longitude  separation  between  two  points,  the  corresponding 
covariance  matrix  is  a  block  matrix  of  Toeplitz  circulant  blocks  (ibid.  pp.  4-5).  Such 
structure  is  exploited  to  enable  a  very  efficient  inversion  of  the  covariance  matrix  as 
demonstrated  by  Colombo  (1981)  in  harmonic  analysis  of  data  on  the  sphere. 


A  Toeplitz  matrix  of  dimension  (N+l)  x  (N+l),  is  a  matrix  Tn  =  [tty],  where  for 
each  element  tty  it  holds  that: 

^kj  =  hc-j  »  k,  j  =  0,  1, ...  ,  N.  (2.3-1) 

Then  the  matrix  may  be  presented  in  the  form: 
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r0 

t-1 

t-2 . 

t.N 

tl 

l0 

t-1 . 

t-(N-l) 

Tn  = 

h 

tl 

(0 . 

t-(N-2) 

*N 

lN-l 

lN-2 

l0 

m 

(2.3-2) 


The  matrix  Tn  is  a  persymmetric  matrix,  that  is,  it  has  a  symmetry  with  respect  to  the 
secondary  diagonal.  There  exists  a  general  class  of  integral  equations  which  reduce  to 
linear  systems  with  Toeplitz  type  matrices  (Tyrtyshnikov,  1980).  The  solutions  of  this 
class  have  the  following  characteristic  properties:  a)  the  kernel  is  invariant  under  some 
transformations,  and  b)  the  region  of  integration  is  obtained  from  parts  of  it  by  applying 
transformations  under  which  the  kernel  is  invariant.  These  invariance  properties 
correspond  to  the  stationarity  property  of  the  stochastic  processes  and  the  associated 
covariance  function.  In  addition  to  the  Toeplitz  structure,  the  covariance  matrices  are 
symmetric  and  positive  definite. 


A  special  case  of  a  Toeplitz  matrix  is  a  circulant  Toeplitz  defined  by  the  condition: 

t-k  =  tN+l-k.  k  =  1, 2, ....,  N.  (2.3-3) 

For  example,  such  matrix  is  the  covariance  matrix  of  data  equally  spaced  over  a  complete 
circle,  when  the  data  belong  to  a  stationary  process.  Also,  the  Toeplitz  circulant  pattern 
appears  for  data  regularly  sampled  completely  over  a  sphere,  as  shown  by  Colombo 
(1979).  The  advantage  in  inversion  of  such  matrices  is  also  explained  by  Eren  (1980),  by 
demonstrating  their  diagonalization  in  the  frequency  domain. 


Of  particular  interest  in  this  study  are  block-Toeplitz  matrices.  Such  structure  is 
completely  equivalent  to  block  matrices  with  Toeplitz  blocks,  which  is  the  one  followed 
by  Colombo  (1979).  It  can  be  proved  that  a  block  Toeplitz  matrix  is  transformed  to  a 
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block  matrix  with  Toeplitz  blocks,  and  vice  versa,  by  reordering  of  the  data 
(Tyrtyshnikov,  1980). 

Consider  the  data  arrangement,  as  shown  in  the  figure  1,  of  equally  spaced 
observations  along  four  parallels  and  three  meridians: 


92 

93 

<p4 


Figure  1.  Data  arrangement  for  the  biock-Toeplitz  structure. 


The  covariance  matrix  of  the  observations  arranged  along  the  meridians  (i.e.  Mi,  M2, 
M3)  may  be  written  in  terms  of  block  matrices  (i.e.  Mj  Mj,  Mi  M2,  ...  etc),  which  are 
teh  covariance  matrices  between  the  observations  on  the  corresponding  meridians  (i.e.  M] 
and  Mi,  Mi  and  M2, ...  etc). 


C  = 


M,Mj 


symmetric 


M,M2 

m2m2 


m,m3 

m2m3 

m3m3 


(2.3-4) 


Then  the  matrices  MjMi,  M2M2,  ...  etc,  according  to  equations  (2.2-34  to  36)  are  as 
follows: 


MjMj  = 


f(<Pi)  f (<Pi><P2)  f(<Pi.<P3)  f((P1,94) 
f((p2)  f(<p2.93j  f(<P2.94) 


symmetnc 


fW  f(<P3.94) 

fW 


(2.3-5) 


The  notation  f  (<pi,  ((>2)  stands  for  "function  of  (pi  and  q>2”  and  so  on.  It  is  easily  seen  that 
the  covariances  among  the  data  within  the  meridians  M2  and  M3  are  the  same  as  for  the 
meridian  Mi,  i.e: 


M1M1  =  M2M2  =  M3M3. 


Now  consider  the  covariances  among  meridians  separated  by  AX. 


M,M2  = 


f  (tp^AX)  f  (tp^AX)  f (tpj.tp^AX)  f (cpj.tp^AX) 

f(i(p2,AX)  f |tp2»<P3> AXj  f(tp2,(p 4, AXj 
symmetric  f(tp3,AX)  f  l^tp^AXj 

f  ((P4*AX) 


(2.3-6) 


(2.3-7) 


It  is  also  obvious  that  the  matrices  reflecting  the  covariances  between  any  two  meridians 
separated  by  AX  are  the  same: 


MiM2  =  M2M3. 


(2.3-8) 


As  already  stated  above  the  matrix  C  in  equation  (2.3-4)  is  a  block  Toeplitz  matrix,  where 
the  blocks  are  general  symmetric  matrices.  Had  the  observations  been  ordered  along  the 
parallels,  instead  of  along  the  meridians,  the  matrix  C  would  have  been  a  general  block 
matrix,  where  the  individual  blocks  are  Toeplitz  matrices.  Note,  that  it  is  not  necessary 
for  the  parallels  to  be  equally  spaced,  but  only  for  the  meridians,  so  that  the  resultant 
matrix  C  is  of  the  form  described  above. 
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It  follows  from  the  Toeplitz  structure  and  the  symmetry  of  the  matrix  C  that  the  top 
row  of  blocks  completely  defines  the  matrix.  In  addition,  the  blocks  are  themselves 
symmetric  matrices,  which  substantially  reduces  the  computations  of  forming  the  matrix, 
as  well  as  the  storage  requirement.  Furthermore,  the  pattern  of  the  inverse  of  Toeplitz 
and  block-Toeplitz  matrices  is  such  that  the  matrix  is  defined  by  one  row,  while  the  rest 
of  the  rows  are  easily  computed  by  recurrence  relationships  (Kutikov,  1967).  The 
algorithm  developed  by  Kutikov  for  the  case  of  block  Toplitz  matrices  is  presented  next. 


Consider  the  block  matrix  of  order  (N+l)  x  (q+1)  defined  as: 


Rn- 


r 


NO 


r 


NN 


(2.3-9) 


The  algorithm  has  been  derived  under  the  following  conditions: 

a)  Rn  is  non-singular, 

b)  Hcs  (for  k,  s  =  0,  1, ...,  N)  are  square  matrices  of  order  (q+1),  such  that  rks  =  rskT, 
and 

n  T 

c)  X  Ckrkscj>0,  (n  =  0,  1, ....  N), 

k.s=0 

i.e.  the  blocks  are  positive  definite,  and  thus  Rn  is  positive  definite. 

In  order  to  describe  the  structure  of  Rn  the  block  matrices  A  and  B  are  introduced, 
defined  from  the  relation: 

Rn  =  BtAB,  (2.3-10) 


which  is  valid  for: 
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—  Si 

O 

and  B  = 

Bio  I 

0 

[  ann 

Bnobni 

I 

where  Ann  and  Bnic  (n,  k  =  0,  1, ... ,  N)  are  square  matrices  of  order  (q+1).  Considering 
block-Toeplitz  matrices,  the  matrix  Rn  of  equation  (2.3-9)  is  written  in  the  form: 


Rn- 


ro  r,  ... 
ri  r0 


r 


N  1  N-l 


(2.3-12) 


Associated  with  the  matrix  Rn  are  the  matrices  Ank  and  Qnn.  introduced  as  the  Bnk  and 
2 

Ann  in  equations  (2.3-1 1).  Then  the  following  recurrence  relations  hold  (ibid): 


®nn  —  Ann  -  E 

n  =  0, 

l,  ....  N 

(2.3-13) 

Bn+1  Jc  =  Bnjc-l  +  Bn+io  Ann.fc  , 

k=  1, 

...,  n  and 

n  =  1, 

...,  N-l 

(2.3-14) 

An+ljc  =  Anjc-1  +  An+1,0  Bn,n-k  » 

k=  1, 

....  n  and 

n  =  1, 

....  N-l 

(2.3-15) 

Bn+1,0  =  *Eni2nn  » 

n  =  0, 

...,  N-l 

(2.3-16) 

2 

An+1,0  =  '  Hn  Ann  , 

n  =  0, 

...,  N-l 

(2.3-17) 

Ann=  X i  BnkE,,.^ 

n  =  0, N 

k=0 

(2.3-18) 

-2  n  T 

^nn=  X  Ankrn  k> 

n  =  0, ...,  N 

k=0 

(2.3-19) 
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En  =  ^  Bnlcrk+1,  n  -  0,  ...,N  1 
k=o 

Hn=  X  ^nk^k+l’  n  =  0,...,N— 1 

k=0 


(2.3-20) 

(2.3-21) 


Now,  let  the  inverse  of  Rn  be  represented  in  terms  of  square  blocks  of  order  (q+1),  then: 


rn  = 


7n 
^ 00 


7n 


7N 

^ON 


7N 

^NN 


(2.3-22) 


Then  they  can  be  computed  from  the  following  equations: 


Bn, 


t  =  0,  1 . N. 


(2.3-23) 


»N 


z;-u-i= Fi-t.s- 


pN 


•.T 


=  zN- 


T  2 


‘N.N-si 


,]« 


NNAN,N-i 


T  2 

+  [bn.s-i]  AnnBn.i-i-  (2.3-24) 

s,  t  =  1,  2, ....  N  and  s  £  t. 


For  the  implementation  of  the  above  algorithm,  the  matrices  B^k,  Ann,  ANk  and  Qnn  are 
computed  using  the  recursive  relationships  (2.3-13)  to  (2.3-21).  Then  equation  (2.3-23) 
is  used  to  compute  the  blocks  of  the  bottom  row  of  Rn  •  Then  all  other  block  rows  are 
computed  by  means  of  equation  (2.3-24).  This  algorithm  is  simplified  for  the  case  of 
symmetric  blocks  Tn  (n  =  0,  1, ....  N).  Thus,  the  formulae  implemented  in  this  work  are 
as  described  above,  but  with 


2  2 
Ann  -  **nn  » 


n  =  0,  1 . N 


(2.3-25) 


Bnk  *  ^nk  » 
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k  =  0,  1, ....  n-1  and 

n  =  1,  2,  ....  N.  (2.3-26) 


In  addition  to  the  Toeplitz  inversion,  a  partitioned  inversion  algorithm  has  also  been 
implemented  in  this  work.  Let 


A  B] 
C  D 


(2.3-27) 


Where  A,  D  and  T  are  square  matrices.  Then  the  inverse  of  T  is  given  (Faddeeva,  1959) 


T  = 


P  Q 
R  S 


where: 


s=(d-ca'1b) 


Q=  -  A  1  BS 
R  =  - SCA  1 


P  =  A  1  -  A  1  BR 


(2.3-28) 


(2.3-29) 


For  a  symmetric  matrix  T  the  above  equations  become: 


A  B 

and  T  1  = 

P  Q 

T  = 

T 

T 

B  D 

Q  S 

(2.3-30) 


where: 
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s  =(d-bta'1b) 

Q=  -  A’1  BS 

-1  -1  T 

P  =  A  -A  BQ 

This  partitioning  is  applied  several  times  for  the  inversion  of  the  covariance  matrix  and  it 
essentially  results  in  a  form  of  stepwise  or  sequential  collocation  (Moritz,  1980,  sec.  19). 
However,  in  that  case,  a  new  term  which  reflects  the  influence  of  new  measurements  is 
added  directly  to  the  predicted  signal  and  the  associated  covariance  matrix. 


CHAPTER  HI 


IMPLEMENTATION  OF  THE  COLLOCATION  PREDICTION  FORMULAE. 


3.1  The  auto-covariance  matrix  of  the  observed  signals 


In  implementing  the  matrix  formulation  (2.2-1 1)  for  the  predicted  potential  harmonic 
coefficients  and  (2.2-14)  for  their  error  covariance  matrix,  one  is  faced  with  the 
theoretical  requirement  of  positive  definiteness  for  the  covariance  matrices.  Additionally, 
the  covariance  matrix  of  the  observed  signal,  Ctt,  must  be  numerically  non-singular  to 
enable  the  computation  of  predicted  signals  for  the  case  of  errorless  observations.  Since 
the  theoretical  and  numerical  properties  of  this  matrix  are  determined  completely  by  the 
covariance  function  and  the  particular  data  sampling  employed,  both  of  these  factors  will 
be  examined  below  with  regard  to  the  positive  definiteness  and  the  singularity  or 
numerical  instability  of  the  matrix  Cu. 


At  first,  the  positive  definiteness  of  the  signal  auto-covariance  matrix  is  examined, 
which  is  equivalent  to  the  positive  definiteness  of  the  covariance  function  in  the  case  of  a 
continuous  signal  (Moritz,  1976).  This  property  is  expressed  by  equation  (2.1-6)  for  the 
covariance  function  K  (P,Q).  Considering  the  form  of  the  fundamental  covariance 
function  implemented  here,  as  given  by  (2.2-20),  this  condition  is  subsequently 
transferred  to  the  degree  variances  Iq,.  By  means  of  the  correspondence  of  the  covariance 
function  and  the  power  spectrum  of  the  disturbing  potential,  the  set  of  all  degree  variances 
kn  constitutes  the  power  spectr  m  (Colombo,  1981,  p.  5).  Then,  a  necessary  and 
sufficient  condition  for  the  covariance  function  to  be  positive  definite  is  the  condition  for 
the  power  spectrum  to  be  positive  (Moritz,  1976,  p.  14),  i.e. 

k„>0,  for  all  n.  (3.1-1) 

The  above  is  also  seen  from  a  different  viewpoint  by  considering  an  equivalent  to  the 
spectral  decomposition  of  the  matrix  Cu,  namely  the  eigenvector  and  eigenvalue 
decomposition  (Colombo,  1979,  p.  14);  the  matrix  of  dimension  q  is  given  in  terms  of 
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the  eigenvalues  X,  and  eigenvectors  m 


i=i 


(3.1-2) 


Then  the  positive  definiteness  of  Ctt  is  assured  for  positive  eigenvalues,  and  the  matrix  is 
invertible  only  for 


Xi>e>0,  fori=l, ....  q  (3.1-3) 

where  e  is  a  sufficiently  large  number  to  avoid  numerical  problems.  Should  the  equality 
be  allowed  in  equations  (3.1-1)  and  (3.1-3),  the  covariance  function  will  be  positive 
semi-definite  and  the  matrix  Cu  will  be  singular,  as  footnoted  by  Moritz  (1980,  p.  93). 
From  the  preceding  arguments  one  can  conclude  that,  when  using  the  series  covariance 
function  expressions,  the  summation  should  theoretically  be  carried  out  to  infinity  in 
order  to  avoid  singularity  of  the  covariance  matrix. 


This  fact  is  also  evident  by  examining  the  matrix  formulation  of  the  collocation 
prediction.  In  particular,  equation  (2.2-13)  where  the  matrix  Cu  is  given  from: 


Cu  =  BKBt. 


(3.1-4) 


As  shown  in  equation  (2.2-48),  K  is  a  diagonal  matrix.  If  the  values  of  kn  are  set  to  zero 
for  n  >  Nmax,  then 

Rank  (K)  =  (Nmax  +  l)2 ,  (3.1-5) 

which  is  the  number  of  coefficients  of  the  spherical  harmonic  expansion  to  degree  Nmax. 
Then,  assuming  that  B  has  full  column  rank, 

Rank  (Cu)  =  Rank  (BKBT)  =  Rank  (K)  =  (Nmax  +  l)2  ,  (3.1  -6) 
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while  the  dimension  of  C«  is  the  number  of  observed  signals.  Finally,  it  may  be  stated 
that  the  disturbing  potential  function  belongs  to  an  infinite  dimensional  space,  thus 
involving  an  infinite  number  of  frequencies,  or  equivalently,  infinite  spherical  harmonics. 
This  fact  is  also  inherent  the  definition  of  the  geodetic  boundary  value  problem  describing 
the  unknown  potential  function  to  be  regular  outside  the  attracting  masses  including 
infinity,  which,  in  turn,  is  by  definition  expandable  to  a  converging  infinite  Taylor  series. 


So  far,  the  covariance  function  analyzed  relates  to  point  data,  as  expressed  by  the 
equations  (2.2-20),  (2.2-26)  and  (2.2-33).  When  the  observed  signals  are  area  mean 
values,  the  corresponding  covariance  function  for  the  gravity  anomaly,  for  example,  is 
given  by  (Sjoberg,  1978,  p.  4). 

covjAgp,  Agq|  =  C(P,  QfdGpd  Cq 

ap°Q  o0  (3.1-7) 

where  ap  and  oq  are  the  areas  on  the  sphere  that  the  values  of  Agp  and  AgQ, 
respectively,  represent.  To  lessen  the  computations  involved  with  the  above  formula,  an 
approximation  is  implemented  in  this  study  as  given  by  equations  (2.2-34)  to  (2.2-36). 
These  formulae  approximate  the  mean  covariance  functions  by  smoothing  the  respective 
point  covariance  functions  by  means  of  the  Pellinen  smoothing  operators  pn  (equation 

2.2-39)).  These  operators,  derived  by  Meissl  (1971,  p.  23),  are  called  smoothing 

operators  in  the  narrow  sense,  since  the  properties  (3n  — >  0  and  Po  =  1  hold,  thus 

n — 

reproducing  the  constant  part  and  damping  the  irregularities  of  the  function  they  are 
applied  to.  This  approximation  of  the  mean  covariance  function  has  been  adopted  by 
Sjoberg  (1979)  for  10°  equal  area  means  successfully,  but  it  is  rejected  by  Colombo 
(1979,  p.  13)  as  too  crude  of  an  approximation,  except  for  the  case  of  a  very  fine  grid. 
The  dependence  of  the  approximation  on  the  area  block  size  is  expected,  since  the 
operators  are  derived  over  a  circular  area  with  angular  radius  \}?o.  which  is  here 
approximated  by  the  block  area  on  the  sphere  (o)  and  the  equivalent  spherical  cap  radius 
q/0,  given  by  (2.2-38). 
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The  possibility  of  utilizing  the  approximation  of  the  mean  covariance  function  with 
the  point  covariance  computed  at  the  appropriate  altitude  as  suggested  by  Tscheming  and 
Rapp  (1974),  has  been  examined  by  Sjoberg  (1978)  and  abandoned  by  him  due  to  the 
poor  results  obtained  in  his  computations. 


Clearly  the  discussion  pertaining  to  the  relation  between  the  point  covariance  function 
and  the  singularity  of  the  covariance  matrix  is  valid  for  the  mean  covariance  functions 
implemented  in  this  work.  They  have  the  same  analytical  form,  but  with  a  smoothed  set 
of  degree  variances  instead  of  kn. 


Regarding  the  actual  computations,  the  truncation  of  the  series  is  necessary,  and 
therefore  a  choice  on  the  maximum  degree  included  in  the  summation  must  be  made. 
Sjoberg  (ibid.)  found  Nmax  =  200  acceptable  for  his  application.  Later  on,  Colombo 
(1981)  analyzed  thoroughly  the  truncation  effect  by  comparing  his  mean  covariance 
function  approximation  (ibid.,  p.  84)  and  formula  (2.2-34)  with  the  rigorous  covariance 
function  in  equation  (3.1-7),  computed  by  numerical  quadratures.  His  experiments 
included  values  of  Nmax  =  180,  300  and  400,  for  5°  and  1°  block  mean  gravity 
anomalies.  The  reported  discrepancies  in  covariance  vary  among  different  cases 
considered.  In  the  case  of  5°  blocks  the  discrepancy  reaches  the  maximum  of  8%  for 
blocks  near  the  pole  and  Nmax  =  180,  whereas  for  1°  blocks  the  difference  is  at  most  5% 
near  the  poles,  for  Nmax  =  180. 

In  this  study  the  values  of  Nmax  =  180  and  360  were  used,  but  the  value  of  Nmax  = 
3000  was  only  used  in  a  single  experiment.  The  sufficient  Nmax  value  was  judged  on 
the  basis  of  the  fi"al  results,  namely  the  recovery  of  the  input  potential  coefficients,  which 
will  be  discussed  in  detail  in  the  next  chapter. 

In  the  following,  figures  2  and  3  show  the  point  covariance  functions  for  Ag  and  N 
respectively,  computed  from  equations  (2.2-26)  and  (2.2-33),  at  1°  intervals,  where  the 
summation  is  carried  out  to  Nmax  =  180.  The  anomaly  degree  variances  cn  have  been 
computed  by  means  of  equation  (2.2-25)  using  the  spherical  harmonic  potential 
coefficient  set  OSU86F.HARMIN.TO360,  described  by  Rapp  and  Cruz  (1986b).  In 


51 


order  to  show  the  covariances  for  the  mean  gravity  anomaly  signal,  figures  4,  5  and  6 
have  been  created  by  plotting  the  point  covariance  function  (as  in  figure  2)  and  the  mean 
covariance  functions,  at  Gfi05  intervals,  for  polar  and  equatorial  blocks,  for  equiangular 
grids  of  10°,  5°  and  3°  respectively.  The  mean  covariances  were  computed  from  equation 
(2.2-34)  with  the  same  parameters  (cn  and  Nmax)  used  in  the  point  covariance 
computation.  It  may  be  seen  that  there  is  a  larger  effect  of  the  smoothing  operators  on  the 
equatorial  block  mean  function,  these  blocks  being  larger  in  area,  while  the  polar  block 
mean  function  rapidly  approaches  the  point  function.  This  is  more  prominent  when 
decreasing  the  grid  size;  in  particular  for  3°  grid  the  difference  between  the  point 
covariance  and  the  polar  block  mean  covariance  is  at  most  17  mgals2.  Similarly,  figures 
7, 8  and  9  show  the  point  and  mean  covariance  functions,  computed  at  (f.05  intervals,  for 
the  undulation  signal.  The  differences  among  the  three  curves  in  each  figure  are  much 
smaller  and  the  effect  of  the  smoothing  is  not  as  prominent,  due  to  the  undulation  being  a 
low  frequency  signal.  In  particular,  even  for  the  5°  grid,  the  polar  block  mean  covariance 
function  almost  coincides  with  the  point  covariance,  with  maximum  difference  of  0.8  m2. 
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Gravity  anomaly  point  covariance  function  computed  from  eq  (2.2-26)  for 
Nmax  =  180  and  degree  variances  implied  b>  .he  OSU86F  field. 
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Figure  4.  Point  and  mean  anomaly  covariance  functions,  for  polar  and  equatorial  blocks 
for  10°  regular  grid,  computed  with  Nmax  =180  and  degree  variances  implied 
by  the  OSU86F  field. 
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Figure  5. 
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Point  and  mean  anomaly  covariance  functions,  for  polar  and  equatorial  blocks 
for  5°  regular  grid,  computed  with  Nmax  =180  and  degree  variances  implied 
by  the  OSU86F  field. 
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Figure  7.  Point  and  mean  undulation  covariance  functions  for  polar  and  equatorial 
blocks  for  10°  regular  grid,  computed  with  Nmax  =  180  and  degree  variances 
implied  by  the  OSU86F  field. 
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Figure  8.  Point  and  mean  undulation  covariance  functions  for  polar  and  equatorial 
blocks  for  5°  regular  grid,  computed  with  Nmax  -  1 80  and  degree  variances 
implied  by  the  OSU86F  field. 
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Figure  9.  Point  and  mean  undulation  covariance  functions  for  polar  and  equatorial 
blocks  for  3°  regular  grid,  computed  with  Nmax  =  180  and  degree  variances 
implied  by  the  OSU86F  field. 
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As  stated  in  the  beginning  of  this  section,  the  data  distribution  is  the  second  factor 
that,  together  with  the  covariance  function,  completely  determines  the  matrix  Cu;  it  is 
discussed  next.  One  related  aspect  is  the  number  of  the  observed  signals,  which  should 
be  sufficient  to  provide  a  solution.  However,  the  signal  should  be  appropriately  sampled, 
since  the  data  spacing  defines  the  frequency  content  of  the  observed  signal  and  therefore 
the  recoverable  information.  In  particular,  for  a  regular  grid  where  Acp  =  AX,  the  global 
sampling  covers  a  system  of  Np  parallels  and  Nm  meridians,  where  Nm  =  2Np,  and  the 
Nyquist  frequency  is  Nq  =  Np.  It  is  known  that  it  is  not  possible  to  estimate  a  complete 
set  of  spherical  harmonic  potential  coefficients  to  degree  and  order  n  >  Nq  (Colombo, 
1981,  p.  11).  When  dealing  with  overdetermined  systems,  as  is  generally  the  case  in 
least-squares  collocation,  repeated  measurements  of  the  same  type,  result  in  identical 
rows  of  the  covariance  matrix  Qt,  and  consequently  in  singularity  (in  errorless 
collocation). 


All  the  above  arguments  point  to  the  fact  that  it  is  neither  necessary  nor  appropriate  to 
include  large  amounts  of  data,  but  instead,  one  should  determine  the  optimal  data 
distribution  which  could  convey  the  maximum  possible  information  to  estimate  the 
gravity  field  up  to  degree  and  order  N  (:  N  <  Nq).  Giacaglia  and  Lundquist  (1972)  have 
derived  an  optimal  grid  of  (N+l)2  sampling  points  on  the  sphere,  by  utilizing  (N+l)2 
independent  sampling  functions.  At  each  point,  only  a  single  sampling  function  has  a 
non-zero  value  and  therefore  the  coefficients  of  this  sampling-function  expansion 
represent  the  gravity  field  at  the  grid  points.  Also,  a  one-to-one  linear  transformation 
between  the  above  expansion  and  the  spherical  harmonic  expansion  to  degree  N  is 
derived  analytically. 


Despite  the  economy  of  such  grids,  the  most  commonly  used  grid  is  the  equiangular 
or  regular  grid,  because  it  is  easily  defined  and  facilitates  the  data  processing.  In  this  type 
of  grid  consisting  of  rows  of  parallels,  each  including  discrete  points,  the  poles  constitute 
singular  points.  Hence,  they  are  eliminated  form  the  grid,  unless  they  are  individually 
handled  as  two  discrete  points,  as  demonstrated  by  Bose  et  al.  (1983).  Still  the  problem 
is  not  completely  overcome,  since  the  polar  areas  cause  numerical  instability  in  the 
covariance  matrix  as  it  tends  to  singularity  for  finer  grids.  This  problem,  also  reported  by 
Colombo  (1979,  p.  13),  is  numerically  shown  in  the  next  two  tables.  The  condition 


numbers  computed  from  the  ratio  of  the  largest  to  the  smallest  eigenvalue  of  the 
covariance  matrices,  are  used  for  a  relative  measure  of  the  effect  of  the  polar  area  data  on 
the  numerical  stability  of  the  matrix. 

The  "limited”  data  sets  are  defined  by  neglecting  data  referring  to  latitudes  larger  than 
±65°  and  ±67.5°  for  10°  and  5°  grids  respectively. 


Table  1.  Condition  numbers  of  covariance  matrices  for  10°  regular  grid. 


Case 

global  Ag 
data  set 

limited  Ag 
data  set 

global  N 
data  set 

limited  N 
data  set 

0 

0 

0 

0 

condition 

number 

0.3x106 

0.9x1 03 

0.2xl0h> 

0.2x106 

Table  2.  Condition  numbers  of  covariance  matrices  for  5°  regular  grid. 


Case 

global  Ag 
data  set 

limited  Ag 
data  set 

global  N 
data  set 

limited  N 
data  set 

number  of 
zero  eigenvalues 

21 

0 

44 

0 

condition 

number 

N/A 

0.2X104 

N/A 

0.4x106 

It  is  seen  that  the  condition  of  the  matrices  generally  improves  when  the  polar  area 
data  are  excluded  from  the  10°  global  data  set  as  indicated  from  the  decrease  of  the 
condition  numbers  by  a  factor  of  103  to  104.  In  the  cases  concerning  the  5°  grid  the 
apparent  singularity  for  the  global  data  sets  is  eliminated  when  considering  the 
corresponding  limited  data  sets. 


Instability  and  singularity  are  also  caused  by  decreasing  the  grid  size,  whereas  the 
global  covariance  function  cannot  distinguish  numerically  between  any  two  adjacent  data 
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points,  thus  producing  two  almost  identical  rows  of  the  covariance  matrix.  Within  the 
same  reasoning  the  undulation  auto-covariance  matrices  are  ill-conditioned  when  the 
corresponding  gravity  anomaly  ones  are  not.  Tables  1  and  2  show  smaller  condition 
numbers  for  the  Ag  covariance  matrices  than  the  corresponding  N  ones.  This  is  due  to 
the  undulation  covariance  function  being  even  less  discriminating,  i.e.  of  low  frequency 
content.  Evidence  to  this  end  may  be  also  found  in  the  figures  presented  in  this  section. 
Specifically  figures  4,  5  and  6  show  a  large  decrease  of  the  covariance  values  for 
spherical  distances  y  :  0°  £  y  £  2°,  while  the  covariances  in  figures  7,  8  and  9  drop 
slightly  within  the  same  interval.  In  addition,  figures  8  and  9  show  that  the  mean 
covariance  for  polar  blocks  coincides  with  the  point  covariance  for  5°  and  3°  regular  grids 
respectively. 


3.2  On  the  regularization  of  the  covariance  matrix. 

Fundamental  definitions  and  concepts  will  be  exposed  below,  aiming  to  provide  a 
concise  interpretation  of  the  experiments  that  will  be  presented  in  the  next  chapter. 

When  solving  a  problem  in  practice,  the  meaning  of  a  "solution"  must  be  defined  as 
well  as  the  features  of  the  computational  algorithm  that  can  be  used  to  obtain  a  solution 

Consider  the  equation 

A  z(s)  =  u  (x)  (3.2-1) 


where  A  is  a  known  continuous  operator,  z(s)  is  the  unknown  function  in  a  space  Si  and 
u  (x)  is  a  known  function  in  a  space  S2-  Suppose  that  for  some  u  =  uj  (x)  the  function 
zi(s)  is  a  solution  of  equation  (3.2-1), 


A  z\  (s)  s  ui  (x). 


(3.2-2) 
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Generally,  values  of  the  function  uj  (x)  are  measured,  and  consequently  only  an 
approximation  u  (x)  of  ui  (x)  is  known.  The  difference  between  u  (x)  and  ui  (x)  is 
measured  in  the  metric  of  the  space  Si  as  the  square  root  of  an  inner  product  psi  (u,,  uj) 
of  the  difference  function  [u  (x)  -  ui  (x)].  An  approximate  solution  of  equation  (3.2-1) 
may  then  be  found,  where  the  error  of  the  approximation  is  measured  in  the  metric  of 
space  S2  by  taking  the  square  root  of  an  inner  product  ps2  (zi,  zj),  of  the  difference 
function  [z  (s)  -  z\  (s)].  For  certain  operators  A  it  can  be  proved  (Tikhonov  and  Arsenin, 
1979,  p.  4)  that  for  an  arbitrarily  small  difference  between  ui  (x)  and  U2  (x),  the 
difference  between  the  corresponding  solutions  zi  (s)  and  z2  (s)  can  be  arbitrary.  In  such 
cases  the  solution  is  not  stable  under  small  changes  in  the  initial  data  u  (x),  which  brings 
up  the  distinction  between  two  groups  of  problems:  well-posed  and  ill-posed  problems. 
It  is  of  interest  here  to  note  that  the  introduction  of  this  concept  was  made  in  the  attempt  to 
develop  appropriate  boundary  conditions  for  differential  equations,  for  example  the 
Dirichlet  and  analogous  conditions  for  elliptical  equations. 

The  precise  mathematical  definition  of  a  well-posed  problem  is  given  below  (ibid., 
p.  7).  Consider  equation  (3.2-1),  where  for  every  element  u  e  S\  with  metric  psi  (ui, 
u2)  there  exists  a  unique  solution 


z  =  R  (u) 


(3.2-3) 


in  the  space  S2  with  metric  ps2  (z\,  z2).  The  metric  is  determined  by  the  formulation  of 
the  problem.  The  algorithm  for  determining  the  solution  z  is  said  to  be  stable  on  the 
spaces  (Si,  S2)  if,  for  every  positive  number  e,  there  exists  a  positive  number  5  (e)  such 
that  the  inequality 

Psi  (ui,u2)£5(e)  (3.2-4) 

implies 


Ps2  (zl>  z2)  -  e  < 


(3.2-5) 
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where  z\  =  R  (ui)  and  Z2  =  R  (112)-  The  problem  of  finding  the  solution  z  in  the  space  S2 
from  the  initial  data  in  the  space  Si  is  said  to  be  well-posed  on  the  pair  of  metric  spaces 
(Si,  S2)  if  the  following  three  conditions  are  satisfied: 

(1)  for  every  element  u  e  Si  there  exists  a  solution  z  e  S2; 

(2)  the  solution  is  unique; 

(3)  the  problem  is  stable  on  the  spaces  (Si,  S2). 

Problems  that  do  not  satisfy  the  above  conditions  are  called  ill-posed.  In  concept,  the 
conditions  (1)  and  (2)  verify  that  the  problem  at  hand  is  determinable  from  the 
mathematical  point  of  view.  Condition  (3)  indicates  whether  the  problem  is  numerically 
determinable  from  the  approximate  data. 


At  this  point  it  is  clearly  understood  that  the  definition  of  an  approximate  solution  of 
an  ill-posed  problem  is  vague.  Still,  this  is  not  a  sufficient  reason  to  completely  avoid 
any  such  problem,  since  there  is  a  large  class  of  problems  for  which  workable  solutions 
are  necessary.  For  example,  the  inverse  operator  A*1  of  the  continuous  operator  A  in 
equation  (3.2-1)  is  not  generally  continuous  on  Si  and,  in  that  case,  the  solution  will  not 
be  stable  under  small  changes  in  u.  Procedures  may  be  developed  for  finding  possible 
solutions  by  utilizing  additional  quantitative  information,  thus  leading  to  a  quasi-solution. 
On  the  other  hand,  when  qualitative  information  is  available,  a  procedure  called 
regularization  can  be  used  to  construct  stable  approximate  solutions. 


The  essential  concept  of  the  regularization  method  is  that  of  the  regularizing  operator 
(ibid.,  p.  47).  An  operator  R  (u,  a)  depending  on  a  parameter  a  is  called  a  regularizing 
operator  for  equation  (3.2- 1 )  if 

(1)  there  exists  a  positive  number  8],  such  that  R  (u,  a)  is  defined  for  every  a  >  0  and 
every  u  e  Si  for  which 


psl  (u,  uj)  <  8  <  81 


(3.2-6) 


with  uj  being  the  exact  value  of  u,  and 
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(2)  there  exists  a  function  a  =  a  (8)  such  that,  for  every  e  >  0,  there  exists  a  number 
5<e)<5i  for  which  the  inequality 

psi  (uj,  us)  <  5  (e) ,  where  us  6  Si  (3.2-7) 


implies 


Ps2  (ZT>  Zo)  ^  e.  where 

(3.2-8) 

Za  =  R  (us,  a  (5)). 

(3.2-9) 

Thus,  a  solution  Z&,  called  a  regularized  solution,  can  be  obtained  from  R  (us,  a  (5)), 
where  a  is  a  numerical  parameter  called  the  regularization  parameter,  by  definition 
consistent  with  the  accuracy  5  of  the  initial  data. 

Although  there  is  no  uniqueness  assumption  for  the  regularizing  operator,  it  is  critical 
to  devise  methods  for  constructing  regularizing  operators.  Usually  the  variational  method 
is  implemented  (ibid.,  p.  50)  where  a  certain  functional,  called  the  smoothing  functional 
is  minimized.  For  example: 

a  2 

M  [z,  u]=  pS)|Az,  u)  +  a  Q(z) ,  (3  2-10) 

with  Q  (z)  a  stabilizing  functional.  Another  useful  method  makes  use  of  the  spectrum  of 
the  operator  A. 


The  method  outlined  in  principle  above,  known  in  literature  as  Tikhonov 
regularization,  was  implemented  by  Rummel,  et  al.  (1979).  Considering  equation  (3.2-1) 
where  u  and  z  belong  in  real  Hilbert  spaces  and  A  is  a  continuous  linear  operator,  a 
minimum  norm  condition  is  invoked.  In  particular: 
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a  >0 


(3.2-11) 


where  uo  is  the  observation  vector,  D  is  the  k  priori  covariance  matrix  for  the  noise  and 
Czz  the  a  priori  covariance  matrix  of  z.  Then,  the  solution  is  obtained  from 


.-t 


aVa  +  oC 


4' 


atd- 


u0 


(3.2-12) 


or  equivalently  from 


ze=C2zA 


(A  CaA 


-l 

+  a  Dl  u 


(3.2-13) 


The  only  difference  between  the  above  equation  and  the  least-squares  collocation  formula 
(2.2-11)  is  the  regularization  parameter  a,  although  in  practice  it  is  often  implied  by 
scaling  the  actual  error  covariance  matrix  D. 


Equation  (3.2-13)  is  recommended  by  Colombo  (1979,  p.  15).  and  it  is  also 
implemented  in  this  study,  when  the  grid  size  is  small  so  as  to  give  rise  to  unstable 
covariance  matrices.  The  instability  is  manifested  both  in  the  elements  of  the  inverse 
matrix  and  in  the  estimated  signal  vector  which  shows  large  discrepancy  from  the  true 
solution.  It  is  important  that  caution  should  be  exercised  in  selecting  the  regularization 
parameter  a,  in  a  way  that  will  stabilize  the  solution  but  also  minimize  the  loss  of 
resolution  in  the  initial  data,  which  is  expected  due  to  the  regularization  process.  For  this 
reason  it  is  beneficial  to  interpret  the  parameter  a  as  a  factor  varying  the  weight  with 
which  the  covariance  information  Czz  influences  the  solution.  The  essence  of  a  may  be 
also  seen  in  the  frequency  domain  by  analyzing  the  eigenvalues  of  the  data  covariance 
matrix  Cuu  =  A  Czz  AT.  Again,  the  fact  that  the  parameter  a  alters  the  spectrum  of  the 
data  covariance  shows  that  it  should  be  chosen  compatible  with  the  magnitude  of  the  error 
in  the  initial  data. 
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3.3  Computational  algorithm  and  software  development 


Since  it  has  been  a  considerable  part  of  this  eifort  to  develop  workable  software  for 
the  boundary  value  problem  solutions  studied  here,  it  is  reasonable  to  present  the 
sequence  of  required  computations,  together  with  the  basic  characteristics  of  the 
algorithms,  the  storage  schemes  and  some  computer  execution  times. 


Most  of  the  programs  were  developed  and  tested  on  the  IBM  3081-D  mainframe 
computer  with  32  megabytes  of  memory,  available  at  The  Instruction  and  Research 
Computer  Center,  at  The  Ohio  State  University.  The  detailed  computations  were  done  on 
a  CRAY  X-MP/48  Supercomputer  available  at  the  Pittsburgh  Supercomputing  Center 
(PSC),  at  Carnegie  Melon  University  and  the  University  of  Pittsburgh. 


Although  no  special  effort  was  made  to  improve  the  source  code  so  as  to  take  full 
advantage  of  the  vector  processing  of  the  CRAY  X-MP/48,  worthwhile  savings  in 
execution  time  were  still  realized.  It  seems  very  realistic  to  expect  that  substantial 
improvement  will  be  achieved  with  certain  software  modifications  to  enable  large  scale 
solutions  and  the  determination  of  high  degree  spherical  harmonic  expansions  for  the 
disturbing  potential.  A  crucial  factor  to  this  is  the  continuing  upgrading  of  the  CRAY 
supercomputers.  The  CRAY  X-MP/48  supercomputer  at  PSC  is  currently  (January, 
1989)  being  replaced  by  the  CRAY  Y-MP/832,  while  a  CRAY-3  is  scheduled  to  arrive  in 
October,  1990.  A  comparison  among  the  three  mentioned  supercomputers  with  regard  to 
their  features  is  given  in  table  3. 
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Table  3.  Supercomputer  upgrading  at  PSC. 


Hardware  specifications 

X-MP/48 

Y-MP/832 

CRAY-3 

Number  of  processing  units 

4 

8 

16 

Memory  in  Megawords 

8 

32 

512 

Clock  cycle  in  nsecs 

9.5 

6 

2 

Floor  space  in  square  feet 

112 

98 

<16 

Performance  in  MFLOPS* 

840 

2,700 

16,000 

*  MFLOPS  "millions  of  floating-point  operations  per  second",  is  the  overall  measure  of 
large-scale  computing  performance.  It  indicates  the  number  of  floating  point  additions, 
multiplications,  etc,  which  the  system  can  do  in  one  second. 

Next,  the  sequence  of  computations  is  described  briefly.  Figure  10  is  offered  for 
further  illustration  of  the  whole  algorithm  in  a  block  diagram  where  each  block  represents 
a  separate  computer  program. 
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Figure  10.  Sequence  of  computations 


In  the  beginning  the  covariance  matrix  is  formed  in  three  steps  using  as  input  a 
spherical  harmonic  potential  coefficient  set  and  the  data  distribution.  The  latter  consists  of 
the  number  of  data  groups,  and  for  each  group,  the  grid  specification  and  location, 
namely:  the  number  of  parallels,  number  of  meridians,  grid  spacing,  latitude  and 
longitude  of  the  utmost  northwest  point,  and  data  type.  During  the  first  step  the  auto¬ 
covariance  matrix  of  the  largest  data  set  is  formed.  This  matrix  is  a  block-Toeplitz  matrix, 
as  described  in  section  2.3;  therefore,  only  the  top  block-row  elements  are  computed,  and 
for  each  block  of  dimension  Np  (Np  is  the  number  of  parallels  in  the  group)  only 
Np(Np+l)/2  elements  are  formed.  A  total  number  of  Nm  (Nm  is  the  number  of 
meridians)  symmetric  blocks  are  formed  and  stored  in  a  vector  array  on  disk.  The 
elements  actually  computed  are  shown  in  the  shaded  area  of  the  diagram  of  figure  1 1 . 


Figure  1 1 .  Sufficient  elements  for  the  auto-covariance  matrix  in  a  single  data  group. 

In  the  second  step,  the  types  of  computations  performed  in  the  first  step  are  done  in  a 
sequential  mode,  and  an  auto-covariance  matrix  is  formed  for  each  data  group.  Then  all 
are  stored  in  one  vector  on  disk.  In  the  last  step,  the  cross-covariances  between  any  two 
data  groups  are  computed.  These  covariances  are  also  formed  into  matrices  which  are 
rectangular  and  therefore  do  not  have  the  Toeplitz  pattern.  However,  some  pattern  exists, 
that  may  be  called  "pseudo-Toeplitz",  since  it  possesses  a  similarity  to  the  Toeplitz  one. 
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This  pattern,  illustrated  in  figure  12,  allows  for  the  matrix  to  be  fully  defined  by  its  first 
block-row  and  first  block-column,  with  each  block  in  position  (i,  j)  being  repeated  at  the 
positions  (i+1,  j+1),  (i+2,  j+2)  etc. 


Np 


— Np2- 


4 


NPlxNmi 


t 


Figure  12.  Sufficient  elements  for  the  data  cross-covariance  matrices. 

Such  structure  can  offer  substantial  savings  in  computer  time,  but  it  has  not  been 
exploited  in  this  work.  After  computing  the  covariances  for  all  possible  combinations  of 
data  groups,  they  are  stored  in  a  vector  on  disk. 

In  the  next  step  of  the  sequence,  the  inversion  of  the  covariance  matrix  takes  place. 
This  is  again  done  in  two  steps,  as  indicated  in  the  diagram  of  figure  10.  First,  the 
covariance  matrix  of  the  largest  group  is  inverted  by  the  Toeplitz  algorithm.  The 
equations  for  this  algorithm  have  been  presented  in  detail  in  section  2.3.  The  input 
consists  of  the  vector  containing  the  auto-covariance  matrix  elements  (shown  in  figure  11) 
and  the  output  is  obtained  in  blocks  (of  dimension  Np  x  Np)  starting  form  the  bottom 
block-row.  When  the  core  memory  is  sufficient  the  elements  of  the  blocks  are  stored  in 
the  matrix  in  symmetric  storage  mode.  For  solutions  requiring  larger  memory  the 
individual  blocks  are  stored  on  disk  as  they  are  computed  or  they  are  multiplied  directly 
with  the  data  vector.  Execution  times  for  forming  the  covariance  matrices  for  gbbal  data 
coverage  and  various  grid  sizes  arc  presented  in  table  4,  together  with  the  r  >pective 
Toeplitz  inversion  times,  for  both  IBM  and  CRAY  compute  runs  when  possible. 
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Table  4.  Time  consumed  in  forming  and  inverting  Toeplitz  covariance  matrices. 
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Next,  the  partitioned  inversion  is  performed  in  a  sequential  mode,  as  illustrated  in 
figure  13.  Each  time  another  group  of  data  is  added,  the  most  recently  computed  inverse 
is  used  in  symmetric  storage  (e.g.  matrix  00)  together  with  the  auto-covariance  matrix  of 
the  next  data  group  (e.g.  matrix  11)  and  the  corresponding  cross-covariances  (e.g.  01). 
Then,  the  augmented  inverse  is  computed  according  to  the  algorithm  given  in  section  2.3 
and  stored  is  symmetric  storage  where  it  can  be  used  as  input  when  a  new  group  is 
considered. 


Figure  13.  The  configuration  of  the  partitioned  inversion. 

The  simulated  data  is  generated  using  a  spherical  harmonic  potential  coefficient  set 
and  the  integrated  associated  Legendre  functions  residing  on  disk.  The  Fortran  programs 
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F428AV1  and  F431AV2  were  used  from  Prof.  R.H.  Rapp's  program  library,  Dept,  of 
Geodetic  Science  and  Surveying.  Program  F428AV1  was  used  for  the  associated 
Legendre  function  computation  utilizing  Paul's  algorithm  (Paul,  1978).  A  modification 
of  program  F431AV2  was  used  to  derive  mean  gravity  anomalies  and  mean  geoid 
undulations  on  a  sphere.  This  program  was  developed  by  Colombo  and  it  is  described  in 
detail  by  Colombo  (1981).  The  Fast  Fourier  Transform  algorithm  implemented  in  this 
program  generates  most  efficiently  the  data  values  along  parallels,  creating  a  global  set  of 
values  referring  to  the  center  of  the  blocks  as  far  as  the  computation  of  the  spherical 
distance,  \j/,  is  concerned.  A  short  program  was  written  to  retrieve  the  appropriate  data 
values  given  a  grid  configuration  and  a  particular  geographic  location  by  defining  the 
utmost  northwest  point.  The  following  program  was  developed  to  collect  all  data  values 
in  a  vector,  add  random  noise,  and  multiply  it  with  the  inverse  of  the  covariance  matrix 
thus  obtaining,  what  is  called  here,  the  first  solution  vector. 


After  a  slight  modification,  program  F428AV1  was  used  to  compute  the  integrated 
associated  Legendre  functions  in  groups  in  accordance  with  the  order  of  the  data  in  the 
observation  vector.  These  data  were  stored  on  disk  and  used  in  the  program  that  was 
designed  to  compute  the  covariance  matrix  between  the  estimated  coefficients  and  the 
data,  on  an  element-by-element  basis.  These  elements  were  multiplied  by  the  appropriate 
elements  of  the  first  solution  vector  and  were  summed  into  two  vector  arrays;  one  for  the 
Cnm  and  one  for  the  Snm  coefficients,  thus  producing  the  estimated  coefficients. 
Substantial  time  savings  were  made  possible  by  using  recursive  formulae  for  the 
computation  of  cos  mX  and  sin  mX  required  in  the  equations  (2.2-58)  to  (2.2-61).  In 
particular,  the  formulae: 


sin  qoc  =  2  sin  ((q- 1 )  a)  cos  a  -  sin  ((q-2)  a) 
cos  qa  =  2  cos  j(q-l)  a)  cos  a  -  cos  M  a) 


(3.3-1) 

(3.3-2) 


were  used,  where  only  three  previously  computed  sines  and  three  cosines  were  saved 
thus  completely  avoiding  any  f,irther  trigonometric  function  calculations.  Another 
version  of  this  program  takes  the  inverse  of  the  data  covariance  matrix  as  input  and 
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computes  the  estimated  coefficients  and  their  estimated  accuracy,  using  one  row  of  the 
cross-covariance  matrix  at  a  time.  In  this  approach  the  entire  inverse  matrix  is  needed  in 
core,  thus  increasing  the  memory  requirements,  while  the  computations  are  overall  very 
time  consuming. 

Finally,  statistics  for  the  recovery  of  the  input  potential  coefficients  are  computed 
using  program  FI 59  from  Prof.  R.H.  Rapp's  library.  Also,  additional  statistics  and 
differences  are  computed  as  discussed  in  detail  in  the  next  chapter. 

The  software  used  in  the  computational  scheme  described  in  this  section  was  written 
by  the  author,  unless  credited  otherwise.  In  particular,  the  source  code  for  forming  the 
various  types  of  auto-  and  cross-covariance  matrices,  the  block-Toeplitz  and  sequential 
inversion  routines  and  the  source  for  the  calculation  of  the  potential  coefficients  and  their 
accuracy  constitute  the  majority  of  the  utilized  software  and  were  developed  and  tested  by 
the  author. 


CHAPTER  IV 


NUMERICAL  RESULTS  AND  ANALYSIS 


4.1  General  experiment  strategy 


Chapters  I  and  II  focus  on  conveying  the  role  of  the  overdetermined  problem  within 
the  general  framework  of  the  geodetic  boundary  value  problems,  and  in  presenting  the 
justification  for  a  collocation  estimation  of  the  spherical  harmonic  parameterization  of  the 
disturbing  potential.  The  formulation  of  the  methodology  is  given  through  the  equations 
of  chapter  n,  while  theoretical  and  practical  aspects  of  the  implementation  are  discussed  in 
chapter  HI. 

This  chapter  is  included  for  the  analysis  of  the  experiments  performed  with  the 
purpose  of  evaluating  the  methodology  in  terms  of  both  accuracy  and  overall 
applicability.  The  experiments  were  designed  to  utilize  synthetic  data  thereby  providing  a 
definite  control  for  the  evaluation  of  the  results  and,  at  the  same  time,  simplifying  the  data 
processing. 


Starting  with  a  given  spherical  harmonic  potential  coefficient  set,  in  particular 
OSU86F.HARMIN.TO360  (Rapp  and  Cruz,  1986b),  point  and  mean  gravity  anomalies 
and  geoid  undulations  are  computed.  Also,  the  anomaly  degree  variances  implied  by  the 
same  field  are  used  in  the  computation  of  the  anomaly  and  the  undulation  auto-covariance 
matrices,  as  well  as  the  cross-covariances  between  the  potential  coefficients  and  the 
anomaly  and  undulation  signals.  The  specific  details  of  the  actual  calculations  are  given 
in  section  4.2.  A  set  of  potential  coefficients  is  estimated  using  collocation  and  the 
simulated  anomaly  and  undulation  data  and  the  necessary  covariances.  The  estimated 
coefficients  are  compared  against  the  input  coefficients,  thus  determining  their  recovery. 
A  number  of  different  quantities  are  then  computed  in  order  to  measure  the  recovery  of 
the  coefficients.  These  quantities  are  the  following: 

( 1 )  The  correlation  between  the  two  coefficient  sets  per  degree  n,  as  computed  from  the 
equation: 
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,  2<n<Nmax 


Rn  = 


Qn 


2  '2 
°n°n 


1/2 


(4.M) 


where 


q»=  X 

m=0 


r  r  +  s  s 

v-nm  '-'nm  J  run  J 


•J 

nm/  * 


2 


O 


n 


=  X(cL+sL) 


and 


m=0 


(4.1-2) 


(4.1-3) 


(4.1-4) 


The  reference  coefficients  are  denoted  by  Cnm  and  Snm,  while  the  estimated  ones  are 
denoted  by  Cnm  and  S^m 

(2)  The  average  correlation  between  the  two  coefficient  sets  is  also  computed  by 
averaging  the  correlations  from  (4.1-1),  i.e: 


R  = 


Nma* 


X  *. 


i=2 

(Nmax  -  1) 


(4.1-5) 


(3)  The  percentage  difference  (PN)  for  each  degree  (n)  is  computed  from  the  equation: 
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PNn  = 


X(acL+as^4 


m=0 
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m=0 
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w 


x  100  , 


(4.1-6) 


where  ihe  difference  between  the  input  and  the  recovered  coefficients  is: 


ACjur,—  Cnm  -  Cnm  and  AS  njn  —  Snm-  S  ^ . 


(4.1-7) 


(4)  The  percentage  difference  is  calculated  for  the  complete  coefficient  set  as  an  average 
given  by: 


PS  = 


Nmax 

I  PN, 

i=2 

Nmax  -  1 


(4.1-8) 


(5)  The  recovery  of  the  input  field  is  also  measured  in  terms  of  the  difference  in  the 
gravimetric  quantities.  In  particular;  the  root  mean  square  difference  in  undulation  per 
degree  is  computed  as: 


UNAn  =  R 


x(acL+as]J 

m=0 


(4.1-9) 


and  cummulatively  to  degree  Nmax: 


A  =  R 


Nmax  n  /  2 

I  I|ACnm+AS 


n=2  m=0 


(4.1-10) 
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where  R  is  the  earth  mean  radius. 

(6)  Similarly  the  root  mean  square  difference  in  gravity  anomaly  is  computed  from: 


AACn  =  y(n-l) 


£(acL+as 

m=0 


and  cummulatively: 


TG  =  y 


-.W 


Nmax  ,2  n  /  2  2  \ 

I  (»-•)  SlACnm+ASj 

n=2  m=0 


(4.1-11) 


(4.1-12) 


(7)  The  recovery  of  the  input  coefficients  is  also  evaluated  with  respect  to  the  estimated 
accuracy  (a)  of  each  coefficient  by  calculating  the  quantities: 


DNr=Sn^Cnjn  ^ 


(4.1-13) 


DNS 


nm 


(4.1-14) 


(8)  Finally,  an  average  value  of  the  above  statistics  is  computed  per  degree  as 


RMS  = 


t  (DNCL 

m=0 


+  dns: 


Ncoef 


(4.1-15) 


where  Ncoef  is  the  number  of  coefficients  in  the  degree  n. 
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Within  the  general  construction  of  the  experiments  as  defined  above,  there  are  several 
inter-related  factors  that  influence  the  problem  and  the  final  recovery  of  the  input 
coefficients.  Since  the  most  critical  objective  of  this  study  is  the  combination  of  different 
data  types,  all  aspects  have  been  analyzed  with  regard  to  the  two  data  types:  gravity 
anomaly  and  undulation.  The  specific  aspects  analyzed  and  presented  in  the  following 
sections  include:  a)  the  representation  of  the  data  as  point  or  mean  values,  b)  the 
approximation  of  the  appropriate  covariance  functions  and  the  truncation  in  the  series 
formulae,  c)  the  error  of  the  data  in  connection  with  the  regularization  procedure  and  d) 
the  size  of  the  regular  grid  and  the  recoverable  frequencies  of  the  field- 


4.2  Specifications  for  the  data  simulation 


This  section  specifies  the  details  of  the  input  in  the  actual  computations  that  were 
outlined  generally  in  the  previous  section. 

Boundary  values  of  mean  gravity  anomalies  and  geoid  heights  were  computed  with 
program  F431AV2  (Colombo,  1981),  which  in  principle  implements  the  equations  (2.2- 
54)  and  (2.2-55).  Modifications  were  made  to  compute  the  data  on  the  surface  of  a 
sphere  of  radius  Rg  using  the  equations: 
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Snm  f  Pnmfcose)sinmXda 

tr 


(4.2-2) 


where  o  is  given  from  equation  (2.2-56)  and 


GM 


(4.2-3) 


Since  the  OSU86F  spherical  harmonic  potential  coefficients  (Rapp  and  Cruz,  1986b) 
were  used  here,  the  associated  scale  factor  with  this  particular  expansion  is  also  used, 
which  is  R  =  ae  =  6378136  m.  In  addition  the  GRS80  constants  were  used,  specifically. 


f  =  1/298.2572 

GM  =  0.3986005  x  1015  mVs2  and 
to  =  7.292115  x  10'5  rad/sec. 

The  normal  field  implied  by  the  GRS80  constants  is  subtracted  from  the  input 
coefficients  by  subtracting  the  zonals  of  degrees  2,  4  and  6.  Then  a  global  set  of  mean 
gravity  anomalies  or  undulations  is  computed.  The  maximum  degree  of  the  expansion 
(Nmax)  and  the  olock  size  are  given  as  input,  although  the  latter  is  limited  to  exact 
dividers  of  tc,  so  that  the  global  number  of  latitude  belts  is  even.  Under  this  scheme  there 
are  no  blocks  crossing  the  equator,  and  the  data  for  two  latitudinal  parallel  rows  that  are 
symmetric  with  respect  to  the  equator  are  calculated  simultaneously,  by  taking  advantage 
of  the  symmetric  and  anti-symmetric  properties  of  the  integrated  Legendre  functions.  The 
data  calculated  for  this  work  refer  to  blocks  of  10°,  5°,  3°  and  2°  with  maximum  degree  of 
expansion  18,  36,  60  and  90  respectively. 


After  the  data  were  generated  and  stored  on  disk  ordered  along  the  parallels,  it  was 
retrieved  in  a  whole  or  partly  and  reordered  along  the  meridians,  in  order  to  obtain  the 
block-Toeplitz  pattern  implemented  in  the  covariance  matrices.  At  the  same  time  the 
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option  of  adding  random  error  is  made  available,  by  utilizing  a  Gaussian  random  deviate 
generator  routine  from  the  IMSL  software  package.  First,  the  subroutine  generates  a  set 
of  uniform  pseudo-random  numbers  in  the  exclusive  range  (0,  1)  given  a  seed  value. 
Then,  these  numbers  are  transformed  to  normal  (0,  1)  deviates  using  the  inverse  normal 
probability  distribution  function.  Finally,  normal  (0,  c2)  deviates  are  computed  by 
scaling  the  generator  output  by  the  standard  deviation  a. 

Consistent  with  the  data  simulation  is  the  computation  of  the  auto-covariance  and 
cross-covariance  matrices  by  means  of  the  equations  given  in  section  2.2.  The  OSU86F 
potential  coefficients  were  used  in  equation  (2.2-25)  to  compute  the  anomaly  degree 
variances  cn  after  removing  the  normal  zonals  of  degree  2,  4  and  6  computed  from  the 
GRS80  constants.  The  coefficient  scaling  factor,  R,  was  used  and  the  value  of  the  ratio 
s,  in  equation  (2.2-26),  was  taken  as  s  =  0.999617.  Then  the  radius  Re  is  defined  as 

R2 

Rr  =  ~ —  —  6379357.8  m  , 

s  (4.2-4) 

which  is  several  kilometers  larger  than  the  earth  mean  radius. 

The  cn  values  as  computed  refer  to  a  sphere  of  radius  R,  but  the  equations  used  f  •*. 
the  computation  of  the  data  covariances  (eg.  equation  (2.2-34))  refer  the  cn  values,  t  d 
therefore  the  covariances,  to  a  sphere  of  radius  Re  by  means  of  the  term  (R2/R£)n+2  as 
explained  by  Rapp  (1988). 

In  case  that  the  radius  R  in  equation  (2.2-34)  is  different  than  the  coefficient  scaling 
factor,  the  anomaly  degree  variances  computed  from  these  coefficients  should  be 
transformed  accordingly  (Moritz,  1980,  p.  181).  This  is  generally  done  to  improve  the 
way  the  series  expression  converges  by  choosing  the  radius  of  an  imbedded  sphere  Re- 
Then  the  cn  values  are  to  be  multiplied  by  the  term  (R/Rb)2i'+4,  which  is  an  option  built  in 
the  programs  for  the  covariance  calculations.  However,  due  to  the  use  of  simulated  data 
only,  the  radius  Re  is  adjusted  instead. 
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4.3  The  effect  of  the  covariance  functions  on  the  potential  coefficient 
recovery 

The  first  issue  to  be  addressed  is  whether  it  is  preferable  to  use  point  or  mean  values 
as  boundary  data.  Sjoberg  (1978)  did  not  examine  this  in  his  work,  because  the  available 
data  base  at  the  time  consisted  of  actual  area  means.  On  the  other  hand,  Colombo  derived 
the  formulae  to  utilize  either  point  or  mean  data  for  harmonic  analysis  of  the  disturbing 
potential  on  the  sphere.  He  underscored  the  complexity  of  this  issue  but  did  not  provide  a 
definite  answer  (Colombo,  1979,  p.  22).  Some  intuitive  argument  was  given  in  support 
of  representing  the  data  as  area  means;  in  this  case  the  aliasing  is  expected  to  decrease, 
since  the  averaging  has  smoothed  out  the  higher  frequencies  and  only  slightly  modified 
the  lower  ones.  No  further  insight  was  made  available,  since  in  the  simulation  studies 
reported  by  Colombo  (1981,  p.  33)  only  area  mean  data  were  considered  "because  area 
means  are  preferred  for  collating  information,  particularly  on  a  global  basis,  at  present". 


Since  no  conclusive  information  is  available  from  previous  research,  it  was 
considered  reasonable  to  undertake  a  first  experiment  with  point  boundary  data.  For  the 
specific  experiment  point  gravity  anomaly  data  were  simulated  at  a  10°  regular  grid  on  the 
sphere  Rg,  using  Nmax  =  18  for  the  degree  of  expansion.  The  auto-covariance  matrix 
was  computed  by  truncating  the  summation  of  the  point  covariance  function  at  180. 
Also,  the  values  of  18  and  36  for  maximum  degree  of  truncation  resulted  in  singular 
matrices  as  seen  from  the  extremely  large  elements  in  the  calculated  inverse  and  the  zero 
eigenvalues.  The  reason  fo-  f  singularity  is  explained  in  section  3.1.  These 
computations  ended  in  absolute  failure  to  recover  the  input  potential  coefficients.  The 
difference  between  the  input  and  estimated  coefficients  was  about  100%  for  most 
degrees,  except  for  the  second  degree  which  was  30%.  These  differences  translate  to 
RMS  undulation  difference  of  40.4  m  and  RMS  anomaly  difference  of  17.8  mgal. 


The  legitimate  question  arises  to  explain  the  above  outcome,  assuming  there  are  no 
software  errors.  In  order  to  test  the  possibility  of  software  errors,  the  computer  program 
written  by  Sjoberg  (1978),  included  in  Prof.  Rapp's  program  library  as  program  F373, 
was  used.  Program  F373  estimates  by  means  of  collocation  the  potential  coefficients  to 
degree  and  order  18  and  their  associated  accuracies.  A  global  set  of  mean  gravity 
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anomalies  is  given  as  input,  which  represent  10°  equal  area  blocks  in  the  system  defined 
by  Hajela  (1975).  The  auto-covariance  matrix  is  derived  by  interpolating  the  tabulated 
covariance  function,  in  particular  Pellinen’s  approximation,  given  by  equation  (2.2-34), 
with  Nmax  =  200.  First,  the  10°  equal  area  means  were  substituted  by  point  values  on  a 
10°  equiangular  grid.  Then,  the  tabulated  mean  auto-covariance  function  was  replaced  by 
the  point  auto-covariance  function,  and  finally,  the  source  code  was  changed  to  compute 
the  cross-covariances  between  the  estimated  potential  coefficients  and  the  point  data. 
After  each  one  of  these  steps,  a  set  of  coefficients  was  estimated.  It  is  clear  that  such 
intermediate  solutions  are  theoretically  incorrect;  nevertheless,  they  indicate  the  sensitivity 
of  the  results  with  respect  to  these  three  components:  the  data,  the  auto-covariance  and 
the  cross-covariance  matrix.  It  was  evident  that  the  change  in  the  cross-covariance  matrix 
drastically  altered  the  estimated  coefficients  to  the  point  that  no  recovery  was  obtained, 
thus  indicating  that  it  has  the  strongest  influence  in  the  transition  from  point  to  mean  data. 

Ai>  a  next  step,  the  point  data  experiments  were  abandoned  and  mean  data 
experiments  were  performed.  The  first  tests  were  set  up  in  correspondence  to  the  point 
data  test,  thus  area  means  of  10°  equiangular  blocks  of  gravity  anomaly  and  undulation 
with  maximum  degree  of  expansion  18  were  used.  The  auto  and  cross-covariance 
matrices  were  computed  accordingly  as  explained  in  section  4.2,  with  maximum  degree 
of  truncation  in  the  auto-covariance  function  of  180.  The  results  of  these  two 
experiments  are  tabulated  in  tables  5  and  6  for  the  Ag  and  N  data  respectively.  In  the 
same  tables  the  results  of  two  other  experiments  are  shown,  where  the  mean  auto¬ 
covariance  matrices  have  been  substituted  by  the  point  ones.  Although  this  configuration 
is  not  recommended  due  to  the  inconsistencies  involved,  they  are  presented  here  in 
support  of  two  arguments;  first,  the  strong  influence  of  the  cross-covariance  matrix  when 
convening  from  point  to  mean  data,  as  observed  with  tests  using  Sjoberg’s  software,  and 
second,  the  distinguishable  behavior  between  the  anomaly  and  undulation  data  and  their 
covariance  functions  with  regard  to  the  recovery  of  the  input  coefficients.  In  this  case  the 
point  covariance  may  be  viewed  as  another  approximation  of  the  mean  covariance 
function  implemented  here. 


Note  that  in  the  tables  presenting  the  recovery,  the  term  "overall"  pertains  to  the 
average  of  the  percentage  difference  (from  equation  (4.1-8))  and  the  RMS  value  for  the 
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undulation  difference  (from  equation  (4.1-10))  and  the  anomaly  difference  (from  equation 
(4.1-12)).  Also,  the  recovery  of  the  degree  17  is  given,  since  the  Nyquist  frequency, 
Nq=18,  of  the  experiments  cannot  be  resolved  completely,  as  explained  in  section  3.1,  p. 
60.  This  fact  is  manifested  by  the  non-recovery  of  the  Nq  frequency  zonal,  i.e.  Cjs  o. 
while  all  other  order  coefficients  within  this  degree  are  recovered  with  some  accuracy. 


Table  5.  Recovery  for  10°  Ag  data,  with  mean  cross-covariance  matrices  and  point  and 
mean  auto-covariance  matrices  with  Nmax  =  180 


Percent 


ESISBUSII 


mean 


36.82 

82.77 

93.49 

92.51 


DX3:] 


2.20 

1.71 

1.75 

3.76 
18.07 
41.87 
36.25 


Undulation  difference  (m)  Anomaly  difference  (mgals) 


mean 


0.416 

0.165 

0.130 

0.215 

0.205 

0.286 

0.254 


2.007 

2.102 

0.938 

0.639 

0.648 


Table  6.  Recovery  for  10°  N  data,  with  mean  cross-covariance  matrices  and  point  and 
mean  auto-covariance  matrices,  with  Nmax  =  180. 


Percent  difference  I  Undulation 


mean 


EiSi4ttua-iuml 


Anomaly  difference  (mgals). 


12 

17 

18 


overall 


45.60 
16.01  24.44 

26.82  86.10 
18.28  17.78 

33.82  9.93 

77.08  45.47 

69.22  29.25 


7.74 


0.485 


10.175 
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It  is  seen  in  table  5,  that  the  Ag  set-up  gives  a  recovery  of  the  input  field,  which  is  30%  to 
40%  better  for  the  low  degrees  (i.e.  for  n  =  2  to  6),  than  for  the  higher  degrees  (i.e.  n=16 
to  18).  As  expected,  in  this  case  the  point  auto-covariance  approximation  is  not  a 
desirable  one,  due  to  resulting  60%  average  difference.  On  the  other  hand,  the  undulation 
data  recovery  is  overall  very  poor,  especially  for  the  low  degrees,  while  the  degree  3  is 
non-recoverable  at  all.  The  degrees  between  10  and  14  are  recovered  somewhat  better 
than  with  the  Ag  data.  Noticable  improvement  is  seen  for  the  case  of  the  point  auto- 
covariance  function,  especially  for  the  degree  3.  This  is  in  agreement  with  the  discussion 
in  section  3.1,  identifying  the  instability  problem  arising  from  the  low  frequency 
undulation  covariance  function. 

The  other  important  issue  addressed  in  this  section  is  the  one  of  the  appropriate 
degree  of  truncation  in  the  auto-covariance  functions,  which  is  discussed  in  principle  in 
section  3.1.  Theoretically  the  summation  should  be  carried  out  to  infinity,  with  the  higher 
frequencies  damped,  according  to  the  block  size.  Colombo’s  tests  (1981,  pp.  85-90) 
have  already  been  summarized  in  section  3.1.  He  reported  that  the  approximation  used 
here  (equation  (2.2-34))  is  inadequate  for  blocks  near  the  poles,  as  compared  to  a  more 
rigorous  covariance  computation  by  means  of  numerical  quadratures.  Such  behavior  of 
the  covariance  function  would  not  be  anticipated,  since  it  only  depends  on  the  block  size 
and  generally  improves  for  finer  grids,  as  is  the  case  near  the  poles.  Also,  the 
covariances  computed  for  this  work  do  not  show  such  discrepancies.  Tables  7  and  8 
present  these  covariances  between  5°  blocks  for  equatorial  and  polar  rows  respectively. 
These  covariance  values  show  similar  agreement  between  the  Pellinen  approximation  and 
the  numerical  quadrature  computation  (columns  1  and  3)  for  both  equatorial  and  polar 
rows,  which  indicates  that  Colombo's  computations  for  the  polar  rows  with  100% 
discrepancy  (column  2,  Table  8)  are  not  reliable.  Note  that  the  cn  values  used  by 
Colombo  were  derived  from  the  "2L"  model  of  Jekeli  (1978),  which  probably  accounts 
for  the  differences  between  Colombo's  and  present  computations  for  the  equatorial  rows 
(columns  2  and  3,  Table  7). 
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Table  7.  Mean  auto-covariance  comparison  of  5°  Ag  data,  reported  by  Colombo  (1981) 
and  present  calculations,  Nmax  =  180,  latitudinal  row  between  0°  and  5°  N. 


Block  No. 

N  umerical/Colombo 

Pellinen/Colombo 

Pellinen/Present 

1 

1  253.$ 

251.8 

2l9.3  H 

2 

149.5 

152.2 

3 

93.9 

93.9 

104.6 

4 

57.1 

57.1 

61.4 

5 

31.7 

31.7 

29.6 

6 

13.9 

13.9 

7.6 

12 

-18.1 

-18.1 

-14.1 

24 

9.1 

9.1 

10.4 

36 

-13.7 

-13.6 

-21.3 

Table  8.  Mean  auto-covariance  comparison  of  5°  Ag  data,  reported  by  Colombo  (1981) 
and  present  calculations,  Nmax  =  180,  latitudinal  row  between  80°  and  85°  N. 


F/VS1  RTStl  r*lH 

Pellinen/Present 

i 

^TT2 

8353 

361.4 

2 

318.3 

800.1 

334.4 

3 

229.2 

709.0 

274.3 

4 

196.4 

592.2 

220.9 

36 

58.0 

149.4 

61.1 

In  this  study  the  overall  recovery  was  used  as  the  criterion  to  judge  the  sufficient 
degree,  Nmax,  of  the  truncation  in  the  summation  of  the  covariance  function  (equations 
(2.2-34)  to  (2.2-36)).  Tests  for  Ag  and  N  were  made  for  Nmax  =  360  and  a  test  for 
Nmax  =  3000  was  performed  for  N  data  only,  since  the  N  type  of  tests  displayed 
instability  with  regard  to  the  covariance  functions  used.  The  input  potential  coefficient  set 
being  complete  to  degree  360,  the  adoption  of  cn  model  was  necessary  to  compute  the 
remaining  cn  values  to  degree  3000.  In  particular,  the  Tscheming-Rapp  model  (1974) 
was  used,  defined  from: 


_  A|n-1) 
C"  (n-2)(n+B) 


n  >  2 


(4.3-1) 
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where  A  =  425.28  and  B  =  24.  The  recovery  obtained  for  the  Ag  data  is  identical  with 
the  one  showed  in  table  5,  thus  proving  the  summation  of  the  covariance  function  to 
Nmax  =  180  satisfactory.  Table  9  contains  the  results  for  the  N  tests,  where  the  symbol 
#  indicates  that  the  value  is  the  same  as  the  one  for  Nmax  =  360  in  the  preceeding 
column. 


Table  9.  Recovery  for  10°  N  data,  mean  covariance  matrices  and  truncation  at  Nmax=360 
and  3000. 


Degree 

Percent  difference 

Undulation  dii 

Ml  W  i  t«l  k'I«  1  ft  H  1 1 1  iTl  W  ■ 

360 

■ 

in 

■ax«M 

MfjSKOm 

X 

# 

rm 

# 

KB 

X 

# 

20.78 

# 

MM 

24.96 

# 

2.409 

# 

1.11 

# 

IHf 

81.14 

81.13 

6.230 

6.229 

3.83 

# 

Kl 

17.86 

# 

1.020 

# 

0.78 

# 

12 

9.95 

# 

0.113 

# 

0.19 

# 

17 

45.50 

# 

0.311 

# 

0.77 

# 

18 

29.25 

# 

0.205 

# 

0.54 

# 

overall 

_ ! 

48.81 

# 

68.818 

11.41 

# 

By  comparing  table  6  with  table  9  it  may  be  seen  that  there  is  a  small  improvement  going 
from  Nmax  =  180  to  360,  but  there  is  no  difference  for  Nmax  =  3000,  which  increases 
tremendously  the  computation  time.  In  all  cases  the  recovery  remains  very  poor,  thus 
making  the  small  difference  completely  trivial. 

4.4  Data  error  consideration  and  regularization  effects. 

The  procedure  for  incorporating  error  in  the  simulated  data  is  presented  in  section 
4.2.  According  to  this  algorithm  the  value  of  the  variance,  c2,  of  the  white  noise  (0,  a2) 
included  in  the  data  must  be  assigned.  The  appropriate  c  is  found  by  means  of  the 
accuracy  estimates  of  the  input  potential  coefficients.  As  explained  in  section  1.2,  the 
estimated  accuracies  of  the  coefficient  set  OSU86D,  described  by  Rapp  and  Cruz 
(1986a),  may  be  used  for  the  OSU86F  field  to  degree  250.  Using  program 
F184.JAN26,  the  error  degree  variances  in  gravity  anomaly,  m£  (Ag),  and  undulation  m„ 
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(N),  are  computed  from  the  estimated  errors  mc^  and  mSnm  of  the  coefficients  Cnfn  and 
Snm  respectively: 

mJ(Ag)  =  Y2{n-l)2(^j  £  (m^-f-  m2J, 

'  E/  m=0  (4.4-1) 

mn(N)=REZ(mL+m*J- 

(4.4-2) 

Then,  the  cumulative  estimated  error  to  degree  Nmax  for  these  two  quantities  is  given 
from: 


Nmax 


X  (As) 

n=2 


and 


oN=± 


Nmax  .  ,  , 

X  mn(N) 

n=2 


1/2 


(4.4-3) 


(4.4-4) 


The  following  table  shows  the  values  of  both  anomaly  and  undulation  standard  deviations 
for  Nmax  =  18,  36,  and  90,  which  correspond  to  the  expansions  estimated  in  this 
work. 


Table  10.  Anomaly  and  undulation  standard  deviations  to  degree  Nmax,  derived  from 
the  estimated  accuracy  of  the  OSU86D. 


HUB 

18 

- 33 - 

— 65 

OAe  (mgals) 

1.15 

2.24 

3.20 

4.61 

on  (meters) 

0.62 

0.79 

0.85 

0.90 
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When  error  is  included  in  the  simulated  observations,  the  associated  variance  is  taken  into 
account  in  the  error  covariance  matrix  designated  by  Cnn  or  D  in  equations  (2.2-1 1)  to 
(2.2-14).  In  order  to  maintain  the  block-Toeplitz  structure  of  the  covariance  matrix  for 
each  data  group,  the  related  D  matrix  is  assumed  diagonal  with  the  additional  assumption 
of  the  same  error  among  observations  within  the  same  parallel.  For  the  computations 
described  here,  the  same  variances  o/g  and  oft  were  assigned  to  all  gravity  anomaly  and 
undulation  observations  respectively.  The  individual  experiments  where  data  error  was 
included  pertain  to  18  and  36  degree  expansion  estimations  and  they  are  described  below. 


The  first  test  utilized  point  anomaly  data,  despite  the  unacceptable  results  obtained  to 
this  point  (Section  4.3).  The  incentive  for  this  was  the  fundamental  difference  in 
truncating  the  degree  of  summation  in  the  point  and  the  mean  covariance  functions. 
When  calculating  the  mean  function,  with  Nmax  =  180,  the  Pellinen  operators  smooth  out 
the  higher  frequencies  according  to  the  data  block  size.  However,  when  calculating  the 
point  function  with  the  same  Nmax  (i.e.  180)  for  10°  point  data  (with  frequency  content 
to  degree  18),  additional  frequencies  in  the  range  18  to  180  are  introduced  in  the  data 
covariances,  while  not  present  in  the  data.  By  including  error  in  the  data,  a  positive 
definite  error  covariance  matrix,  D,  is  added  to  the  data  covariance  matrix,  thus  making 
the  matrix  inversion  possible  even  for  low  degree  (eg.  18)  truncation  of  the  covariance 
function.  Using  10°  point  Ag  data  with  maximum  degree  of  expansion  18,  the  same  (i.e. 
18)  maximum  degree  in  the  covariance  function  summation  and  a=±1.15  mgals,  the 
results  showed  no  recovery,  as  before. 

The  point  data  configuration  was  abandoned  once  again,  and  the  next  tests  were 
made  with  mean  gravity  anomaly  and  mean  undulation  data  computed  from  an  expansion 
to  degree  18.  The  auto-covariance  matrices  were  computed  with  input  parameters  Nmax 
=  180,  $Ag  =  ±  1.15  mgals  and  3n  =  ±  0.62m  accordingly.  Although  the  adopted 
values  of  ^Ag  and  3n»  as  shown  in  Table  10,  pertain  to  the  error  associated  with  point 
quantities,  they  are  used  in  the  present  analysis  as  rough  estimates  of  the  data  error  in  the 
mean  anomaly  and  mean  undulation  observations.  The  results  obtained  for  the  potential 
coefficient  recovery  are  presented  in  Tables  1 1  and  12. 
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Tabic  11.  Recovery  for  10°  Ag  data  with  error,  mean  covariance  matrices  and 
Nmax=180. 


W*nirth'f‘La 

Undul.  diff.  (m) 

2 

iiMi  1 1  m 

195 

iMnrtT3M 

3 

0.9990 

4.90 

0.926 

0.28 

4 

0.9987 

5.29 

0.511 

0.24 

5 

0.9991 

4.37 

0.323 

0.20 

6 

0.9994 

3.52 

0.201 

0.15 

12 

0.9821 

24.98 

0.283 

0.48 

17 

0.9286 

46.39 

0.317 

0.78 

18 

0.9391 

43.38 

0.304 

0.79 

i—  i7.K 

“  1 

The  effect  of  data  error  in  this  case  may  be  seen  by  comparing  the  results  of  tables  5  and 
11.  The  error  in  the  estimated  coefficients  of  second  degree  increases  from  0.58%  to 
5.95%  when  data  error  is  included.  It  is  seen  that  the  error  increase  is  smaller  for  higher 
degree  coefficients,  while  the  coefficients  near  the  Nyquist  frequency  are  recovcered  with 
an  error  of  about  40%.  The  average  error  increase  to  degree  18  is  about  30%  for  using 
data  with  error.  This  figure  corresponds  to  increase  of  the  RMS  undulation  difference  by 
102%  and  increase  of  the  RMS  anomaly  difference  by  26%. 


Table  12.  Recovery  for  10°  N  data  with  error,  mean  covariance  matrices  and  Nmax=180. 


Degree 

P55HEESM 

I 

1 

mmsb 

2 

2432 

3 

53.79 

X 

3.13 

4 

WfflTfl  -  ' 

19.65 

1.896 

5 

29.66 

2.196 

1.35 

6 

0.9828 

18.48 

1.055 

0.81 

12 

0.9738 

22.80 

0.258 

0.44 

17 

0.7074 

72.47 

0.496 

1.22 

18 

0.7907 

63.04 

0.442 

1.15 

P  33.94 

11.659 

4.69 

A  drastic  improvement  is  apparent  between  the  results  of  Table  1 2  and  the  corresponding 
errorless  case  shown  in  Table  6,  especially  for  degrees  2  and  3,  which  amounts  to  a 
decrease  in  RMS  undulation  difference  by  85%.  This  result  provided  more  evidence  of 
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the  instability  of  the  data  solutions  and  the  regularization  procedure  was  considered  as 
the  most  justifiable  next  step.  The  fundamental  theory  of  the  regularization  is  briefly 
described  in  section  3.2  including  formula  (3.2-13)  which  is  implemented  here.  As 
already  mentioned  the  error  covariance  matrix,  D,  is  assumed  diagonal  in  all  experiments 
presented  here.  Therefore,  the  term  aD  in  equation  (3.2-13)  may  be  written  as  ao2I, 
where  I  is  the  unit  matrix.  In  all  undulation  error  covariance  matrices  from  this  point  on, 
the  value  aa2=9m2  is  used.  For  the  case  of  10°  N  data,  the  corresponding  error 
suggested  in  Table  10  is  5n  =  ±0.62m,  thus  defining  the  value  of  a  to  be  about  23. 
Tables  13  and  14  contain  the  results  for  the  potential  coefficient  recovery,  obtained  from 
two  experiments  with  10°  N  data  and  regularized  auto-covariance  matrices  by  a  factor 
0=23.  In  the  first  experiment  the  noise  (0,  On2)  is  included  in  the  simulated  data  (table 
13),  while  in  the  second  one  (table  14),  errorless  data  are  used. 

Table  13.  Recovery  for  10°  R  data  with  error,  mean  covariance  matrices  with 
Nmax=180  and  regularization  of  a  =  23. 


The  comparison  of  Table  13  to  Table  12  shows  improvement  in  the  recovery  of  the  low 
degree  coefficients,  but  poorer  recovery  of  the  coefficients  of  the  second  half  of  the 
estimated  spectrum.  Thus,  the  average  percent  difference  changes  from  34%  to  36% 
when  regularization  is  applied.  On  the  other  hand,  the  RMS  undulation  difference  is 
decreased  by  16%  and  the  RMS  anomaly  difference  is  decreased  by  6%. 
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Table  14.  Recovery  for  10°N  errorless  data,  mean  covariance  matrices  with 
Nmax=180  and  variance  ac5N=9m2 


■•r.T7T4nir.T.« 

Anomaly  diff1.  (mgals) 

2 

2L91 

3 

0.9835 

42.87 

8.112 

2.49 

4 

0.9930 

16.15 

1.559 

0.72 

5 

0.9807 

23.85 

1.766 

1.09 

6 

0.9898 

14.29 

0.816 

0.63 

12 

0.9934 

0.368 

0.62 

17 

0.7687 

80.13 

0.548 

1.35 

18 

0.9258 

73.95 

0.518 

1.35 

overall 

H2£Z££9i 

r"““  35.44 

9.517 

4.32 

The  overall  improvement  in  the  recovery  when  neglecting  the  error  in  the  simulated  data, 
while  already  including  the  regularized  error  covariance  matrix  is  generally  small,  of  the 
order  of  3%.  However,  a  comparison  on  a  per-degree  basis  of  the  numbers  shown  in 
Tables  13  and  14  shows  that  the  improvement  of  the  overall  recovery  does  not  mean 
improvement  in  individual  degree  recovery. 


Before  presenting  the  tests  for  various  regularization  parameter  values,  the  tests 
including  data  error  in  36  degree  expansion  estimation  are  presented.  Specifically,  mean 
anomaly  and  mean  undulation  data  referring  to  5°  blocks  were  derived  with  maximum 
degree  36.  Both,  anomaly  and  undulation,  auto-covariance  matrices  were  singular  due  to 
decreased  data  spacing.  To  enable  the  inversion  of  the  covariance  matrix  the  value  of  3 
mgals2  was  added  to  the  diagonal  elements  in  the  case  of  Ag  data,  and  the  value  of  9m2  in 
the  case  of  N  data.  Then,  the  noise  added  to  the  Ag  simulated  observations  is  distributed 
as  (0,  1.732)  and  the  noise  added  to  the  N  simulated  observations  is  destributed  as  (0, 
0.792),  which  implies  a  regularization  factor  of  a  =  14.  In  the  following.  Tables  15  and 
16  show  the  results  of  the  Ag  data  tests  with  data  error  and  without  data  error 
respectively.  However,  the  same  error  covariance  matrix  is  included  in  both  cases. 
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Table  15.  Recovery  for  5°  Ag  data  with  error,  mean  covariance  matrices  with 
Nmax=180  and  variance  a3^g2  =  3  mgals2. 


Table  16.  Recovery  for  5°  Ag  errorless  data,  mean  covariance  matrices  with 
Nmax=180  and  variance  aS^g2  =  3  mgals2. 


■ 

I 

■nunruciTaroi 

E53555S!  gESjOKBEBHB 

2 

032 

031 

3 

1.0000 

0.51 

0.03 

4 

1.0000 

0.43 

■ 

0.02 

5 

1.0000 

0.57 

0.03 

6 

1.0000 

1.09 

0.05 

12 

0.9996 

6.33 

0.072 

0.12 

18 

0.9975 

10.92 

0.076 

0.20 

24 

0.9959 

24.44 

0.095 

0.34 

3o 

0.9434 

47.61 

0.142 

0.77 

Ieeemb 

r~  1(5.96 

Including  data  error  increased  the  error  of  the  estimated  coefficients,  as  is  seen  by 
comparing  the  results  of  Tables  15  and  16.  Especially,  the  recovery  of  the  low  degrees  is 
mostly  affected,  for  example  the  error  of  estimating  the  second  degree  coefficients 
increased  from  0.32%  to  3.05%.  In  the  overall  sense,  the  average  error  to  degree  36 
increased  by  25%,  while  the  RMS  undulation  difference  increased  by  104%  when  error  is 
included,  thus  reflecting  the  large  impact  on  the  low  degree  coefficient  recovery. 
However,  the  increase  in  RMS  anomaly  difference  is  about  18%.  It  may  be  pointed  out 
that  the  effect  of  the  data  error  is  generally  smaller  than  the  same  effect  for  the  10°  Ag 
data,  presented  earlier  in  this  section.  The  corresponding  results  for  5°  N  data  with  and 
without  errror  are  presented  in  Tables  17  and  18  respectively. 
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Tabic  17.  Recovery  for  5°  N  data  with  error,  mean  covariance  matrices  with 
Nmax=180  and  variance  ocon2=9  m2 


BU333 

BBS  52BSS51 

K2BSnElS053!!S 

Undul.  diff.  (m) 

Anomaly  diff.  (mgals) 

2 

0.9847 

1731 

"  3.137  ^ 

3 

0.9912 

24.86 

4.703 

1.45 

4 

12.49 

1.206 

0.56 

5 

0.9934 

11.90 

0.881 

0.54 

6 

0.9956 

9.36 

0.535 

0.41 

12 

11.21 

0.127 

0.22 

18 

0.9895 

30.85 

0.216 

0.57 

24 

0.9175 

71.33 

0.278 

0.98 

36 

0.7475 

89.05 

0.267 

1.43 

[■STOf 

MKBI 

P“  502 

5.i?7 

Table  18.  Recovery  for  5°  "FT errorless  data,  mean  covariance  matrices  with 

Nmax=180  and  variance  cc3n2=9  m2. 


■T*7773MiWi« 

2 

TTimrrBl 

15787 

HKS  ITS ■ 

UAA 

3 

0.9906 

23.55 

1.37 

4 

11.48 

MHRi 

0.51 

5 

BP  (79 

11.69 

0.866 

0.53 

6 

8.60 

0.491 

0.38 

12 

0.9995 

8.44 

0.096 

0.16 

18 

0.9965 

29.21 

0.205 

0.54 

24 

0.9761 

67.92 

0.264 

0.93 

36 

0.8999 

88.34 

0.264 

1.42 

Emma 

43.37 

6.651 

“j 

Table  18  shows  better  recovery  for  errorless  data  as  compared  to  Table  17.  In  particular, 
the  error  of  the  second  degree  coefficients  increases  by  10%,  while  the  average  error  to 
degree  36  increases  by  3%.  Note  that  the  large  influence  on  the  low  degree  recovery 
observed  in  the  Ag  tests  is  not  observed  here,  due  to  the  large  regularization  factors.  The 
RMS  undulation  difference  increased  by  7%  while  the  RMS  anomaly  difference  increased 
by  2%.  These  figures  show  an  error  influence  of  the  same  order  of  magnitude  as  the  10° 
N  tests,  presented  in  Tables  13  and  14. 


Finally  seven  regularized  solutions  for  N  errorless  data  are  presented  for  various 
regularization  factors.  Since  generally  the  solutions  improve  for  the  lower  degrees  but 
deteriorate  for  the  higher  frequencies,  they  are  compared  with  each  other  in  terms  of 
overall  recovery  to  various  degrees,  in  particular  4,  9  and  18.  These  statistics  are  given  in 
Tables  19  and  20,  together  with  the  recovery  measures  for  degree  2.  Note,  that  for  the 
computation  of  the  covariance  matrices  the  covariance  function  was  truncated  at 
Nmax=360.  According  to  the  tests  presented  in  section  4.3,  the  results  would  be 
identical  if  the  value  of  Nmax  =  180  were  used  instead. 


Table  19.  RMS  undulation  differences  (in  meters)  for  10°  N  errorless  data  and  mean 
covariance  matrices  (Nmax  =  360)  from  various  regularized  solutions. 
Variance  =  a  a&. 


variance  (m2) 

n  =  2 

n  =  2  to  4 

n  =  2  to  9 

n  =  2  to  18 

4 

4.1*4 

^  1 1.172 

11.2i2 

9 

3.962 

9.358 

9.634 

9.724 

3.582 

7.741 

7.982 

8.174 

50 

2.951 

6.268 

6.669 

100 

2.344 

4.734 

5.119 

5.776 

300 

1.413 

2.980 

4.375 

5.425 

600 

1.031 

4.982 

6.05 

Table  20.  RMS  anomaly  differences  (in  mgals)  for  10°  N  errorless  data  and  mean 
covariance  matrices  (Nmax  =  360)  from  various  regularized  solutions. 
Variance  =  a  aft. 


variance  (m2) 

n  =  2 

n  =  2  to  4 

n  =  2  to  9 

n  =  2  to  18 

4 

wn 

437 

9 

0.61 

4.30 

20 

IlMHrM 

4.54 

50 

ES9 

■SZM 

2.28 

5.04 

100 

0.36 

1.34 

2.33 

5.58 

300 

0.22 

0.86 

3.36 

6.70 

600 

0.16 

0.73 

4.39 

7.51 

The  results  in  Tables  19  and  20  show  a  continuing  improvement  in  the  recovery  of 
degree  2  for  increasing  variance  values.  In  fact,  the  RMS  undulation  difference  is 
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decreased  by  75%  between  variance  =  a  On  =  4  and  variance  =  a  On  =  600.  However, 
the  recovery  of  the  higher  degrees  of  the  estimated  expansion  deteriorates  when  the 
variances  increases.  This  is  seen  more  clearly  by  observing  the  overall  RMS  anomaly 
difference  which  increases  for  larger  variances.  Such  results  are  to  be  anticipated  in 
applying  a  regularization  procedure,  especially  for  large  regularization  factors  as  in  the 
presented  experiments  where  the  value  of  a  varies  from  10  to  1.5  x  103. 

4.5  Recovery  in  relation  to  the  grid  size. 

The  sets  of  spherical  harmonic  potential  coefficients  estimated  from  various  grid  size 
data  are  analyzed  in  this  section.  All  solutions  are  regularized  except  for  the  one  obtained 
from  10°  mean  anomaly  data  to  degree  18.  The  same  value  was  added  to  the  diagonal  of 
the  auto-covariance  matrix  for  each  data  type,  regardless  of  the  block  size.  In  particular, 
the  value  of  3  mgals2  was  used  to  regularize  the  Ag  solutions  which,  according  to  the 
standard  deviations  given  in  table  10,  would  actually  correspond  to  underscaling  the 
variance  of  the  data.  To  the  contrary,  the  results  shown  in  tables  19  and  20  suggest  large 
regularization  factors  in  the  case  of  N  data.  The  value  of  9m2  was  used  in  all  the 
solutions  with  N  data,  which  corresponds  to  regularization  parameters,  a,  in  the  range 
from  10  to  23. 

The  data,  consisting  of  Ag  and  N  values,  were  computed  on  regular  grids  of  size 
10°,  5°,  3°  and  2°  with  maximum  degrees  of  expansion  of  18,  36,  60  and  90  respectively. 
The  auio- covariance  matrices  were  derived  with  the  mean  covariance  function  summed  to 
Nmax  =  180.  Among  the  recovery  statistics,  the  percent  difference  and  the  correlation 
per  degree  are  shown  in  figures  14  to  17  and  figures  18  to  21  for  all  four  solutions 
derived  from  Ag  data.  Similarly,  the  results  obtained  from  N  data  are  shown  in  figures 
22  to  29. 
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Figure  15.  Percent  difference  per  degree  between  the  input  coefficients  and  the 
coefficients  recovered  from  5°  Ag  errorless  data  to  degree  36. 
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Figure  17.  Percent  difference  per  degree  between  the  input  coefficients  and  the 
coefficients  recovered  from  2°  Ag  errorless  data  to  degree  90. 
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Figure  19.  Correlation  per  degree  between  the  input  coefficients  and  the  coefficients 
recovered  from  5°  Ag  errorless  data  to  degree  36. 
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Figure  21.  Correlation  per  degree  between  the  input  coefficients  and  the  coefficients 
recovered  from  2°  Ag  errorless  data  to  degree  90. 
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Figure  23.  Percent  difference  per  degree  between  the  input  coefficients  and  the 
coefficients  recovered  from  5°  N  errorless  data  to  degree  36. 


Figure  25.  Percent  difference  per  degree  between  the  input  coefficients  and  the 
coefficients  recovered  from  2°  N  errorless  data  to  degree  90. 
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Figure  26.  Correlation  per  degree  between  the  input  coefficients  and  the  coefficients 
recovered  from  10°  N  errorless  data  to  degree  18. 
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Figure  27.  Correlation  per  degree  between  the  input  coefficients  and  the  coefficients 
recovered  from  5°  N  errorless  data  to  degree  36. 
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Figure  29.  Correlation  per  degree  between  the  input  coefficients  and  the  coefficients 
recovered  from  2°  N  errorless  data  to  degree  90. 
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Several  observations  may  be  made  from  examining  figures  14  to  29:  (1)  the  lower 
degrees  in  any  estimated  set  of  coefficients  are  recovered  better  than  the  higher  degrees 
within  the  same  expansion,  (2)  the  recovery  of  the  input  coefficients  improves  with 
decreasing  grid  size  for  both  data  types  and  (3)  the  recovery  obtained  from  the  Ag  data  is 
better  than  the  one  obtained  from  N  data. 


Average  recovery  statistics  are  presented  to  degrees  18,  36,  60  and  90.  These 
include  average  correlation,  average  percent  difference,  RMS  undulation  difference  and 
RMS  anomaly  difference  as  defined  in  section  4.1.  Tables  21  to  24  contain  these 
statistics  for  the  Ag  solutions,  and  tables  25  to  28  the  corresponding  results  for  the  N 
solutions. 

Table  21.  Average  correlation  to  various  degrees,  Nmax,  between  the  input  coefficients 
and  the  coefficients  recovered  from  Ag  errorless  data  given  on  various  grids. 


"\block 

nmc 

Nmax 

10° 

D 

a 

H 

18 

0.9918 

0.9994 

1.0000 

36 

0.9929 

Vfffl 

0.9999 

60 

0.9994 

90 

0.9947 

Table  22.  Average  percent  difference  to  various  degrees,  _Nmax,  between  the  input 
coefficients  and  the  coefficients  recovered  from  Ag  errorless  data  given  on 
various  grids. 
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Table  23.  RMS  undulation  difference  (in  meters)  to  various  degrees,  Nmax,  between  the 
input  coefficients  and  the  coefficients  recovered  from  Ag  errorless  data  on 
various  grids. 


Table  24.  RMS  anomaly  difference  (in  mgals)  to  various  degrees^Nmax,  between  the 
input  coefficients  and  the  coefficients  recovered  from  Ag  errorless  data  on 
various  grids. 


Table  25.  Average  correlation  to  various  degrees,  Nmax,  between  the  input  coefficients 
and  the  coefficients  recovered  from  N  errorless  data  given  on  various  grids. 
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Table  26.  Average  percent  difference  to  various  degrees,  _Nmax,  between  the  input 
coefficients  and  the  coefficients  recovered  from  N  errorless  data  given  on 
various  grids. 


^^block 

Nmax  \ 

10° 

5° 

3° 

2° 

18 

35.44 

14.46 

6.62 

4.18 

36 

43.37 

26.28 

15.47 

60 

49.53 

35.59 

90 

54.20 

Table  27.  RMS  undulation  difference  (in  meters)  to  various  degrees,  Nmax,  between  the 
input  coefficients  and  the  coefficients  recovered  from  N  errorless  data  given 
on  various  grids. 


fill 

10° 

5° 

3° 

2° 

18 

9.517 

5.533 

2.940 

2.450 

36 

5.652 

3.028 

2.486 

60 

3.171 

2.590 

90 

2.683 

Table  28.  RMS  anomaly  difference  (in  mgals)  to  various  degrees,  Nmax,  between  the 
input  coefficients  and  the  coefficients  recovered  from  N  errorless  data  given 
on  various  grids. 


"\block 

N^ize 

Nmax 

10° 

5° 

3° 

2° 

18  1 

4.32 

2.03 

1.06 

0.83 

36 

5.28 

3.35 

2.06 

60 

7.63 

5.77 

90 

9.83 

The  general  trend  indicated  by  the  numbers  in  tables  21  to  28  may  be  pointed  out. 
First,  for  both  data  types,  there  is  an  improvement  in  the  recovery  of  the  potential 
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coefficients  with  denser  data.  For  example,  the  improvement  in  average  percent 
difference  up  to  degree  18  is  87%  and  in  RMS  undulation  difference  is  86%,  when 
decreasing  from  10°  to  3°  blocks  of  Ag  data.  The  corresponding  figures  for  N  data  are 
81%  and  70%.  In  decreasing  from  10°  to  2°  blocks  these  numbers  improve  to  95%  and 
92%  for  the  Ag  data  and,  88%  and  75%  for  the  N  data. 

It  is  also  apparent  that  superior  results  are  obtained  using  Ag  data  than  when  using  N 
data.  From  the  statistics  in  the  preceeding  paragraph  the  improvement  with  decrease  in 
grid  size  is  also  superior  for  the  Ag  data.  Consider  the  percent  difference  between  the 
recovered  and  the  input  coefficients  shown  in  tables  22  and  26  for  Ag  and  N  data 
respectively.  The  differences  between  columns  4  of  these  two  tables  are  3.45,  13.41, 
30.96  and  54.20  for  corresponding  Nmax  values  of  18,  36,  60  and  90.  Thus,  the 
difference  in  the  recovery  statistics  between  the  two  data  types  increases  with  degree, 
especially  for  degrees  larger  than  36,  as  the  two  last  differences  indicate.  Also,  Figures 
28  and  29  show  that  the  percent  differences  for  N  data  are  more  than  40%  for  degrees 
higher  than  36. 

In  addition  to  the  coefficients  their  accuracy  is  estimated  by  means  of  the  equation 
(2.2-14)  for  the  expansions  to  degrees  18  and  36.  Due  to  the  large  size  of  the  auto- 
covariance  matrices  for  the  case  of  3°  and  2°  mean  data  (actually  7,200  and  16,200 
observations  respectively)  as  well  as  due  to  the  time  consuming  computation  required,  the 
accuracy  estimates  of  these  expansions  have  not  been  computed. 

Negative  variances  were  computed  for  certain  coefficients  which  generally  suggest 
that  the  covariances  used  are  not  optimal  in  the  sense  of  representing  the  true  covariances. 
In  particular,  the  variances  estimated  from  10°  Ag  data  for  the  coefficients  C30,  C31  and 
S3 1  are  -1.8  x  10*14,  -6.9  x  lO'1^  and  -6.9  x  10'15  respectively.  These  values  are 
practically  zero,  since  the  internal  representation  of  the  double-precision  numbers  retains 
15  significant  digits.  No  negative  variances  were  calculated  in  estimating  the  36  degree 
field  from  the  5°  Ag  data.  However,  several  negative  variances  were  calculated  from  the 
N  data,  especially  for  the  low  degrees,  which  indicate  inadequacy  of  the  covariance 
function  used  for  the  undulation  data.  Tables  29  and  30  show  the  particular  coefficients 
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with  negative  estimated  variances  as  well  as  the  extremes  of  these  variance  values.  The 
symbol  "#"  means  the  same  variance  with  the  one  of  the  Cnm  coefficient. 


Table  29.  Extreme  values  of  negative  variances  and  all  related  coefficients,  derived  from 
10°  N  errorless  data. 


Order 

Largest  Variance 

Smallest  Variance 

mm 

mm 

Cnm 

Snm 

0  to  2 

1,2 

-0.9x10- 11 

-0.3xl0-12 

-0.3xl013 

-0.3xl0-13 

0  to  3 

1  to  3 

-O.lxlO'10 

-0.3xl0-n 

-0.3xl0-14 

-0.3xl0-14 

0  to  4 

1  to  4 

-0.2xl0-12 

-0.5xl0'13 

-0.3xl0-17 

-0.3xl0'17 

0  to  4 

1  to  4 

-0.9xl0"13 

-0.2xl0'13 

-0.3xl0'15 

-0.3x10- 15 

0  to  5 

1  to  5 

-O.lxlO-13 

-0.6x10- 14 

-O.lxlO15 

-O.lxlO'15 

0  to  6 

1  to  6 

-0.3x10' 14 

-0.3xl0-14 

-0.9xl0'16 

-0.9xl0'16 

8 

lto3,  5to7 

lto3,5to7 

-0.2xl0'15 

-0.2xl0-15 

-0.2xl016 

-0.2xl0'16 

mm 

2, 3,6,7, 8 

2,3,6,7,8 

-0.4xl015 

-0.4x10-15 

-O.lxlO'16 

-O.lxlO'16 

■9 

2,6,7 

2,6,7 

-0.2xl0-15 

-0.2xl015 

-0.2xl0'16 

-0.2xl0'16 

Table  30.  Extreme  values  of  negative  variances  and  all  related  coefficients  derived  from 
5°  N  errorless  data. 


Order 

Largest  Variance 

Smallest  Variance 

Cnm 

99EHIH 

Cnm 

beih 

mm 

nan 

-1.9xlO-13 

# 

-1.6xlO'15 

# 

mm 

IS; 

-3.3xlO'13 

# 

-9.  lx  10- 14 

# 

Hi 

mm 

1.  2 

-2.9xl0"14 

# 

-1.3xlO*15 

# 

E9 

1,  2 

-l.lxlO"14 

# 

-7.2xl0'16 

# 

all 

1,  2 

1,  2 

-4.6x  10"15 

# 

-3.1xl0"16 

# 

mm 

■SB 

1,  2 

-2.3xIO"15 

# 

-1.7xl0'16 

# 

H 

1 

-2.7xl016 

# 

N/A 

N/A 

■9 

111 

1 

-9.4xl0'17 

# 

N/A 

N/A 

The  last  recovery  statistics,  i.e.  the  differences  between  the  recovered  and  the  input 
coefficients,  were  computed  and  normalized  by  the  estimated  standard  deviations  of  the 
particular  coefficients  (see  equations  (4.1-13)  and  (4.1-14)).  Then  the  RMS  value  of  the 


110 

normalized  differences  was  computed  for  each  degree  (from  equation  (4. 1-15))  and  was 
plotted  together  with  the  extremes  of  the  absolute  values  of  the  normalized  differences 
within  the  particular  degree.  These  plots  are  presented  for  the  variances  of  the 
coefficients  derived  from  10°  and  5°  Ag  data  in  figures  30  and  31.  Also,  figures  32  and 
33  show  these  results  derived  from  10°  and  5°  N  data.  Note  that  the  negative  variances 
have  been  set  to  zero  and  the  corresponding  coefficient  differences  have  not  been  included 
in  the  RMS  computation. 


Figure  30.  RMS  and  extreme  absolute  values_of  the  normalized  differences  between  the 
input  and  the  recovered  coefficients  from  10°  Ag  errorless  data. 
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Figure  32.  RMS  and  extreme  absolute  values  of  the  normalized  differences  between  the 
input  and  the  recovered  coefficients  from  10°  N  errorless  data. 
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Figure  33.  RMS  and  extreme  absolute  values  of  the  normalized  differences  between  the 
input  and  the  recovered  coefficients  from  5°  N  errorless  data. 

From  figures  30  and  31  it  is  seen  that  the  RMS  values  are  generally  centered  between 
the  extremes  for  each  degree,  and  also  they  increase  with  increasing  degree.  The  RMS 
values  as  well  as  the  ranges  of  the  normalized  differences  are  decreased  for  the  degrees  up 
to  18  when  using  the  5°  Ag  data  instead  of  10°  Ag  data.  However,  the  other  end  of  the 
spectrum  (i.e.  for  24  <  n  <  36)  shows  increased  values  compared  to  the  differences  for 
12  <  n  <  18  in  figure  30. 


Figures  32  and  33  do  not  indicate  as  clear  as  Figures  30  and  31  the  trend  of 
increasing  value  with  increasing  degree,  although  the  RMS  values  tend  to  stabilize  for  the 
second  half  of  the  predicted  specirum.  Also,  the  differences  referring  to  the  predicted 
coefficients  from  the  N  data  are  generally  larger  than  the  corresponding  ones  predicted 
from  the  Ag  data.  Finally,  it  is  clear  that  the  RMS  values  of  the  normalized  differences 
are  smaller  than  one,  i.e.  the  overall  differences  of  the  predicted  coefficients  from  the 
reference  ones  are  within  the  accuracy  of  the  prediction.  There  are  only  three  exceptions, 
as  it  can  be  seen  in  figure  32,  specifically,  values  between  1.5  and  2.0  for  n  =  5,  6  and  8. 
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Similar  graphs  have  been  produced  for  the  experiments  where  data  noise  was 
included,  which  were  discussed  in  section  4.4.  Since  the  solutions  presented  in  this 
section  are  regularized  (i.e.  with  the  error  covariance  matrix  effectively  added  to  the  signal 
covariance  matrix),  the  estimated  accuracy  of  the  predicted  signals  remains  identical  to  the 
errorless  data  case,  as  it  is  clearly  seen  from  equation  (2.2-14).  However,  the  estimated 
coefficients  themselves  are  affected  by  the  data  error,  thus  changing  the  values  of  the 
normalized  differences.  The  coefficient  error  normalized  by  the  corresponding  estimated 
accuracy  is  shown  per  degree  in  Figures  34  and  35  for  5°  Ag  and  N  noisy  data 
respectively. 


O 


0  8  16  32  40 


HARMONIC  DEGREE 


Figure  34.  RMS  and  extreme  absolute  values  of  the  normalized  differences  between  the 
input  and  the  recovered  coefficients  from  5°  Ag  data  with  noise. 
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Figure  35.  RMS  and  extreme  absolute  values  of  the  normalized  differences  between  the 
input  and  the  recovered  coefficients  from  5°  N  data  with  noise. 

By  comparing  Figure  34  to  Figure  31  it  is  seen  that  only  the  first  half  of  the  predicted 
spectrum  is  affected  by  the  data  noise,  thus  bringing  the  normalized  differences  to  the 
same  order  of  magnitude  for  all  degrees,  about  0.5.  Such  effect  is  not  observed  for  the  N 
data.  The  comparison  of  Figure  35  to  Figure  33  shows  only  slight  changes  in  the  low 
degrees.  For  the  other  degrees  there  are  only  small  changes  in  the  extreme  values,  while 
the  RMS  values  are  unaffected. 

4.6  Solutions  combining  data  types. 

So  far,  ail  the  computations  were  made  in  order  to  analyze  the  behavior  of  the  two 
data  types  individually  in  estimating  spherical  harmonic  coefficients  with  collocation. 
This  section  presents  experiments  made  with  combined  data  types,  namely  mean  gravity 
anomalies  and  mean  undulations. 
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Two  different  tests  were  made  to  estimate  potential  coefficients  to  degree  18  and  36, 
using  mean  data  given  on  10°  and  5°  regular  grids,  respectively.  The  10°  data  were 
derived  from  the  OSU86F  expansion  to  maximum  degree  18,  and  the  5°  data  were 
derived  from  the  same  expansion  to  degree  36.  The  auto-covariance  matrices  were 
computed  by  summing  the  covariance  function  to  Nmax  =  1 80.  The  value  of  9m2  was 
added  to  the  diagonal  elements  of  the  undulation  auto-covariance  matrices  and  the  value  of 
3  mgals3  was  added  to  the  anomaly  covariance  matrix  of  5°  block  data  only,  in  order  to 
perform  experiments  comparable  to  the  ones  presented  in  section  4.5. 


Five  data  groups  were  included  in  each  experiment.  The  first  group  is  a  global  Ag 
data  set,  as  used  in  the  experiments  presented  previously  in  section  4.5.  The  other  four 
data  groups  contain  N  data  and  cover  the  following  geographic  areas: 

(1)  -10°  S  ^  <p  <  40°  N  and  140°  E  £  X  £  240°  E,  which  includes  5  rows  of  parallels  and 

10  zones  of  meridians  on  a  10°  regular  grid, 

(2)  0°  <.  <p  ^  60°  N  and  300°  E  <,  X  <  350°  E,  which  includes  6  rows  of  parallels  and  5 
zones  of  meridians, 

(3)  -70°  S  <,  <p  £  -30°  S  and  0°  <,  \  <  120°  E,  which  includes  4  rows  of  parallels  and  12 
zones  of  meridians,  and 

(4)  -60°  S  £  <p  <  - 10°  S  and  180°  E  <  X.  S  290°  E,  which  includes  5  rows  of  parallels  and 

1 1  zones  of  meridians. 

The  data  location  is  also  shown  on  the  map  of  figure  35. 

When  estimating  potential  coefficients  to  degree  18  from  10°  grid  data,  the  number  of 
observations  contained  in  each  group  above  is  50,  30,  48  and  55,  for  a  total  of  183 
undulation  observations.  There  are  648  anomaly  observations,  which  constitute  78%  of 
the  831  total  observations.  In  the  case  of  5°  grid,  the  number  of  data  is  exactly  doubled 
by  densifying  the  coverage  at  the  same  geographic  locations. 
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Figure  36.  Undulation  data  coverage  for  the  combination  solutions. 

Using  the  above  data  configuration,  the  coefficients  were  estimated  by  employing  the 
sequential  type  of  algorithm  which  was  described  in  section  3.3.  The  statistics  of  the 
recovery  of  the  reference  potential  coefficients  resulting  from  these  tests  are  presented  in 
the  following  figures  and  tables.  Figures  37  and  38  show  the  percent  differences  per 
degree  for  the  recovery  obtained  from  10°  and  5°  data  respectively,  while  figures  39  and 
40  show  the  correlation  per  degree  between  the  input  coefficients  and  the  ones  obtained 
from  these  tests. 
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Figure  38.  Percent  difference  per  degree  between  the  input  coefficients  and  the 
coefficients  recovered  from  5°  Ag  and  N  combined  errorless  data  to  degree  36. 
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Figure  40.  Correlation  per  degree  between  the  input  coefficients  and  the  coefficients 
recovered  from  5°  Ag  and  N  combined  errorless  data  to  degree  36. 
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Figures  37  and  38  show  the  same  general  trend  as  figures  14  and  15  respectively,  which 
reflect  recovery  from  Ag  data  alone.  A  more  detailed  comparison  indicates  that  for 
degrees:  2  £  n  <  6  the  percent  differences  in  the  combination  solution  increase  about  1%, 
and  for  degrees:  1  l^n<15  decrease  by  1%,  while  they  are  the  same  for  the  remaining 
degrees.  This  may  be  interpreted  as  the  influence  of  the  incorporation  of  the  N  data  on 
the  Ag  obtained  solution.  Figures  39  and  40  compare  in  a  similar  way  to  figures  18  and 
19. 

Specific  values  of  the  computed  statistics  for  selected  degrees  are  presented  in  table 
31  for  the  18  degree  solution  and  table  32  for  the  36  degree  solution.  The  overall 
recovery  to  degrees  1 8  and  36  is  given  in  table  33. 


Table  31.  Recovery  statistics  for  10°  Ag  and  N  combined  errorless  data. 


Table  32.  Recovery  statistics  for  5°  Ag  and  N  combined  errorless  data. 


orrelation  1  Percent  difference  I  Undul.  diff.  (m) 


1.0000 

1.0000 

1.0000 

.9999 

.9996 

.9974 

.9435 


0.064 

0.077 

0.142 


Anomaly  diff.  (mgals) 


1 
5 
3 

0.02 

0.06 

0.11 

0.20 

0.76 
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Table  33.  Overall  recovery  to  degree  Nmax  obtained  from  Ag  and  N  combined 
errorless  data  given  on  10°  and  5°  regular  grids. 


Nmax 

Grid 

size 

Average 

Correlation 

Average 
%  diff. 

&MS  undulation 
diff.  (m) 

RMS  anomaly 
diff.  (mgals) 

18 

13.84 

1.47 

18 

5> 

4.27 

36 

5> 

BEEJSBK 

16.73 

2.08 

By  comparing  the  second  row  in  table  33  with  the  corresponding  statistics  in  tables  21  to 
24  it  is  evident  that  all  statistics  are  slightly  better,  except  for  the  RMS  undulation 
difference,  which  is  increased  by  3.5  cm. 

Finally,  the  estimated  accuracy  of  the  predicted  coefficients  was  computed.  Several 
negative  variances  were  found  for  the  case  of  10°  data  solution,  in  particular  for  degrees  2 
to  5  and  various  orders,  as  shown  in  table  34.  A  single  negative  value  was  computed  for 
the  C30  coefficient  estimated  from  5°  data,  which  was  of  the  magnitude  -4.7  x  10*15. 

Table  34.  Extreme  values  of  negative  variances  and  all  related  coefficients  derived  from 
10°  Ag  and  N  combined  errorless  data 


Order 

^EEEESHTil  V'f  HEW 

Smallest  Variance 

Hi 

mm 

wM 

Cnm 

S  nm 

ma 

0,2 

2 

-0.2xl0'13 

-0.5xl0-l4 

-0.6X10-14 

N/A 

HI 

0  to  3 

1  to  3 

-O.lxlO'^ 

-0.5xl0-l4 

-0.1x10-13 

-0.3x1 0-I4 

Hi 

3 

3 

-0.5x10-15 

-0.3x10-15 

N/A 

N/A 

Hi 

4 

4 

-1. OxlO-16 

-0.8x10-16 

N/A 

N/A 

In  addition,  the  estimated  standard  deviations  are  used  to  normalize  the  difference  of 
the  respective  estimated  coefficients  from  the  reference  ones,  excluding  the  coefficients 
given  in  table  34.  The  range  of  the  absolute  values  of  the  normalized  differences  as  well 
as  the  RMS  value  for  each  degree  is  shown  in  figures  41  and  42. 
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Figure  42.  RMS  and  extreme  absolute  values. of  the_normaIized  differences  between  the 
input  and  the  recovered  coefficients  from  5°  Ag  and  N  combined  errorless  data. 
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The  comparison  of  figures  41  and  42  with  figures  30  and  31  shows  larger  RMS 
values  and  larger  extremes  for  the  first  half  of  the  predicted  spectrum  in  the  combination 
solutions,  while  the  second  half  seems  unaffected.  On  the  other  hand  the  comparison  of 
these  figures  with  figures  32  and  33  shows  the  same  general  trend  for  the  low  degrees, 
but  with  the  extreme  values  supressed  when  combining  data  types. 


CHAPTER  V 


CONCLUSIONS 


Numerical  solutions  to  the  overdetermined  geodetic  boundary  value  problem  have 
been  derived  using  the  method  of  least-squares  collocation.  The  boundary  value  problem 
is  defined  by  the  Laplacian  for  the  disturbing  potential  and  the  boundary  conditions  of  the 
first  and  third  (or  mixed)  type  holding  on  overlapping  parts  of  the  spherical  boundary  of 
radius  Rg.  The  solutions  are  expressed  in  terms  of  solid  spherical  harmonics,  thus  sets 
of  coefficients  are  obtained  to  the  degree  defined  by  the  Nyquist  frequency  of  the 
boundary  data.  The  boundary  data  are  mean  values  referring  to  equiangular  blocks  on  the 
sphere  Re  and  they  have  been  simulated  using  the  OSU86F  potential  coefficient  set.  In 
particular,  mean  gravity  anomaly  and  mean  undulation  data  were  computed  on  regular 
grids  of  10°  and  5°  using  expansions  to  degree  18  and  36  respectively.  The  equations 
applied  to  estimate  the  coefficients  and  their  accuracies  are  given  in  section  2.2,  as  well  as 
the  equations  used  to  compute  the  data  covariances  and  the  cross-covariance  between  the 
estimated  coefficients  and  the  data.  Specifically,  for  the  computation  of  the  data 
covariances  equations  (2.2-34)  to  (2.2-36)  are  used,  where  the  anomaly  degree  variances 
are  computed  from  the  OSU86F  coefficients  and  the  summation  is  carried  out  to  Nmax  = 
180.  As  a  consequence  of  the  implemented  isotropic  and  homogeneous  covariance 
functions,  the  auto-covariance  matrix  of  data  given  on  a  regular  grid  (defined  by  a  number 
of  parallels  and  meridians)  is  identified  as  a  block-Toeplitz  matrix.  Substantial  savings  in 
computing  effort  are  realized  in  forming  and  inverting  such  matrices.  An  algorithm  is 
implemented  to  exploit  this  pattern  for  the  largest  data  group,  while  additional  groups  are 
included  in  a  sequential  algorithm  as  described  in  section  3.3.  Since  only  simulated  data 
are  analyzed,  the  error  of  the  solution  is  determined  by  means  of  the  recovery  of  the 
reference  field.  The  errors  of  the  estimated  solutions  to  degree  18  and  36  increase  with 
degree,  ranging  from  0.5%  to  maximum  of  10%  for  the  first  half  of  the  estimated 
spectrum  and  reach  40%  near  the  Nyquist  frequency. 


To  consider  the  influence  of  the  two  data  types,  a  number  of  solutions  were  carried 
out  based  exclusively  on  gravity  anomaly  data  or  undulation  data.  By  comparing  the 
combination  solutions  with  the  corresponding  one-data-type  solutions,  the  former  can  be 
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viewed  as  a  weighted  average  of  the  latter.  This  averaging  is  done  in  an  optimal  way 
under  the  least-squares  collocation  minimum  principle,  and  the  combination  solutions 
ultimately  depend  on  the  relative  amounts  of  data  as  well  as  their  distribution.  In  the  case 
that  data  error  is  included,  it  also  affects  the  combination  in  a  relative  sense:  the  more 
accurate  data  has  greater  influence  on  the  combined  solution.  Similar  effects  happen  due 
to  the  regularization  process.  When  large  regularization  factors  are  used,  the  data  error  is 
effectively  increased  and  the  weight  it  carries  on  the  solution  is  decreased.  The  solutions 
presented  in  this  work  are  obtained  with  a  factor  of  10  for  scaling  the  undulation 
variances.  The  anomaly  variances  are  not  scaled,  thus  the  final  result  is  influenced  more 
by  the  anomaly  data,  as  discussed  in  section  4.6. 


Since  the  result  of  data  combination  may  be  predicted  based  on  solutions  from 
anomaly  data  and  solutions  from  undulation  data,  several  parameters  of  the  system  can  be 
analyzed  with  respect  to  the  individual  data  types.  This  does  not  only  simplify  the 
computations,  but  also  provides  a  better  insight  into  the  behavior  of  the  different  data,  and 
some  guidelines  towards  optimizing  their  combination. 

Given  the  analytical  form  of  a  covariance  function,  the  two  fundamental  properties  of 
positive  definiteness  and  symmetry  are  required  in  order  for  it  to  represent  a  true 
covariance  function,  and  also  in  principle,  a  reproducing  kernel  in  the  Hilbert  space  where 
the  collocation  solution  is  derived.  For  the  functions  implemented  in  this  study  the 
positive  definiteness  is  assured  when  the  power  spectrum,  defined  by  the  degree 
variances  kn,  is  positive  for  all  n.  As  it  is  discussed  in  detail  in  section  3.1,  truncation  of 
the  covariance  function  summation  at  finite  degree  n  results,  theoretically  in  singular 
covariance  matrices.  Numerical  tests,  presented  in  section  4.3,  have  shown  that 
truncation  at  Nmax  =  180  is  sufficient  for  data  given  on  a  10°  regular  grid.  Note  that  this 
condition  holds  for  point  and  mean  covariance  functions.  The  issue  of  using  point  vs 
mean  data  has  been  addressed  in  section  4.3.  However,  meaningful  solutions  were 
derived  when  using  area  mean  values  only.  Point  data  on  a  10°  regular  grid  completely 
failed  to  recover  the  input  coefficients.  Other  than  with  regard  to  the  covariance  matrix 
singularity,  the  value  of  Nmax  was  tested  with  regard  to  the  recovery  of  the  input 
coefficients.  The  results  shown  in  tables  5,  6  and  9  prove  that  the  value  of  Nmax  =  180 
is  sufficient,  since  there  is  no  noticeable  change  for  Nmax  =  360  and  3000. 
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After  studying  corresponding  solutions  obtained  separately  from  global  anomaly  and 
undulation  data  sets,  several  statements  may  be  made  about  the  existing  differences  and 
similarities.  First,  the  undulation  point  covariance  function  used  is  generally  a  low 
frequency  function  with  correlation  length  of  about  22°,  whereas  the  correlation  length  for 
the  anomaly  point  function  is  about  2°,  as  is  seen  in  Figures  3  and  2.  These  correlation 
lengths  increase  for  mean  covariance  functions,  especially  with  decreasing  grid  size  as 
shown  in  figures  4  to  9.  As  a  result  of  the  implemented  global  covariance  function 
characteristics,  ill-conditioning  appears  in  the  auto-covariance  matrices,  which  introduces 
instabilities  in  the  solutions.  This  type  of  effect  occurs  in  both  types  of  covariance 
matrices,  although  more  acutely  in  the  undulation  one,  as  a  result  of  using  data  in  polar 
areas  of  an  equiangular  grid,  where  the  blocksize  decreases  when  approaching  the  poles. 
For  data  given  on  a  regular  grid  of  size  equal  or  smaller  than  5°,  singularity  occurs  in  the 
covariance  matrix  of  both  data  types,  as  is  indicated  by  their  condition  numbers  given  in 
tables  1  and  2. 

The  singularity  can  be  handled  by  adding  the  error  covariance  matrix,  D  =  a2  I,  to 
the  signal  covariance  matrix.  In  this  case  error  may  be  included  in  the  simulated 
observations,  in  the  manner  explained  in  section  4.1,  by  adding  pseudo-random  numbers 
distributed  as  Gaussian  (0,  o2).  The  effect  of  the  data  error  is  tested  for  mean  data 
referring  to  10°  and  5°  equiangular  blocks  and  it  is  discussed  in  section  4.4.  When  using 
10°  Ag  noisy  data,  the  average  error  of  the  estimated  coefficients  to  degree  18  increased 
by  30%,  as  opposed  tc  the  errorless  data.  Especially,  the  error  of  the  lower  degrees  is 
mostly  affected,  as  the  results  of  Table  1 1  indicate.  However,  the  results  obtained  from 
lCP  N  noisy  data  are  not  comparable  (Table  12).  In  this  case,  a  regularization  method  is 
recommended,  where  the  variance  of  each  observation  is  scaled  by  a  factor  a>l. 
Specifically,  the  variance  of  a(5&  =  9m2  was  added  to  the  diagonal  elements  of  the  N 
signal  covariance  matrix.  Then,  the  recovery  of  the  input  coefficients  was  tested  for 
noisy  as  well  as  errorless  simulated  observations.  Tables  13  and  14  show  3%  increase  of 
the  average  error  to  degree  18,  due  to  observation  noise.  The  same  increase  was 
observed  in  the  average  error  to  degree  36,  derived  from  5°  N  noisy  data  (Tables  17  and 
18).  The  value  of  3  mgals2  was  used  for  the  variance  of  5°  Ag  data  and  was  added  to  the 
Ag  signal  covariance  matrix.  When  including  the  equivalent  noise  in  the  observations  the 
average  error  to  degree  36  increased  by  25%,  as  the  results  in  Tables  15  and  16  indicate. 
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Tests  made  using  undulation  signals  on  a  10°  regular  grid,  showed  that  the  error  of  the 
estimated  low  degree  coefficients  decreases  with  increasing  regularization  factors,  as 
presented  in  tables  19  and  20.  Specifically,  the  coefficients  of  degree  2  are  estimated  with 
an  RMS  error  of  6%  for  a  =  1.5  x  103,  while  for  a  =  10,  this  error  is  about  30%. 
However,  it  is  known  that  the  regularization  effectively  decreases  the  data  resolution,  i.e. 
increases  the  data  error,  thus  the  error  in  the  estimated  higher  degrees  increases  rapidly 
with  the  factor  a.  The  errors  corresponding  to  the  two  previous  a  values  are  99%  and 
80%  near  the  Nyquist  frequency.  An  indication  of  the  implemented  undulation 
covariance  function  being  problematic,  i.e.  not  approximating  well  the  true  covariance, 
are  the  unrealistic  error  estimates  obtained  by  the  least-squares  collocation  formula;  in  this 
case,  negative  variances  are  calculated  for  coefficients  of  degrees  2  to  10.  These  values 
are  practically  zero  and  are  given  in  tables  29  and  30.  The  data  densification  by 
decreasing  the  grid  size  has  a  positive  influence  on  the  results.  The  improvement  is  well 
manifested  for  both  data  types,  and  it  is  discussed  in  the  comparative  analysis  of  section 
4.5,  by  analyzing  the  overall  recovery  statistics  given  in  tables  21  to  28.  For  example, 
the  average  percent  difference  to  degree  18  is  reduced  by  95%  for  the  anomaly  solutions 
and  80%  for  the  undulation  solutions,  going  from  10°  to  2°  blocksize.  Further 
improvement  is  expected  when  approaching  the  limit  of  continuous  data.  Finally  the 
recovery  ability  of  the  two  data  types  is  compared.  The  anomaly  data,  of  all  block  sizes, 
produce  coefficients  with  error  less  than  1%  for  the  low  degrees,  approaching  a 
maximum  of  10%  for  coefficients  of  degree  half  the  Nyquist  frequency  and  reaching  40% 
near  degrees  of  the  Nyquist  frequency.  The  error  definitely  increases  with  increasing 
degree.  On  the  other  hand,  the  undulation  data  recover  the  2nd  and  3rd  degree  coefficients 
with  an  error  of  10%  -  20%.  Then  the  error  decreases  to  about  5%  for  coefficients  to 
degree  10  and  afterwards  it  continues  to  increase  reaching  90%  near  the  Nyquist 
frequency.  Specific  numbers  are  presented  and  discussed  in  section  4.5. 


As  expected,  the  solutions  require  considerable  computer  time,  even  in  a  CRAY 
supercomputer.  High  degree  solutions  for  a  global  data  set,  although  time  consuming, 
are  manageable  under  the  Toeplitz  covariance  matrix  scheme.  However,  the  time 
requirements  increase  rapidly  when  the  sequential  algorithm  is  implemented  in  adding 
groups  of  data.  In  particular,  for  the  experiments  performed,  the  time  required  for  the 
estimation  of  the  36  degree  expansion  using  a  global  data  set  is  about  2  minutes,  while  the 
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combination  solution  to  degree  36  with  4  additional  data  groups  requires  about  60 
minutes.  Other  than  the  required  computer  time,  the  need  in  storage  is  another  limiting 
factor.  For  example,  the  90  degree  expansion  involves  16,200  observations  and  the 
associated  covariance  matrix.  There  is  no  need  for  storing  this  matrix,  unless  the  error 
estimates  are  computed.  In  this  case,  the  matrix  can  only  be  stored  on  tape  by  writing 
parts  of  it  at  a  time. 

To  end  this  discussion,  several  recommendations  are  made  for  future  improvement 

(1)  The  undulation  global  covariance  function  must  be  studied;  an  improvement  should 
be  made,  which  will  affect  the  solution  in  a  way  similar  to  the  regularization  procedure, 
thus  improving  the  prediction  of  the  low  degree  coefficients. 

(2)  For  further  improvement  of  the  low  degrees,  a  satellite-derived  field  should  be 
included.  This  can  be  done  easily  in  including  these  potential  coefficients  as  an  additional 
data  group  in  the  sequential  algorithm,  which  is  simplified  considering  that  the  errors  of 
the  satellite  coefficients  are  not  correlated  with  the  other  data  errors. 

(3)  Programming  effort  should  be  dedicated  in  order  to  take  advantage  of  the  pseudo- 
Toeplitz  pattern  (section  3.3)  in  forming  the  cross-covariances  between  any  two  data 
groups,  as  well  as  utilizing  these  matrices  in  the  partitioned  inversion  scheme.  Also,  fully 
vectorized  software  should  be  used  to  take  advantage  of  the  continually  improving 
supercomputer  capabilities.  Considering  that  the  CRAY  Y-MP/832,  currently  in 
operation  at  PSC,  is  at  least  3  times  faster  and  has  4  times  larger  memory  than  the  CRAY 
X-MP/48  used  in  this  work,  combination  solutions  to  degree  90  are  feasible.  Depending 
on  the  improvement  that  can  be  achieved  by  implementing  the  suggestions  made  here  and 
on  the  impact  that  fully  vectorized  software  will  have,  high  degree  solutions  to  degree  180 
may  be  possible.  However,  when  judging  the  required,  time  of  the  methodology  of  this 
study  in  comparison  with  other  methods  utilizing  one  data- type,  the  effort  required  for  the 
data  type  conversion  should  be  considered. 
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